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A METHOD OF CHARACTERISTICS FOR STEADY 
THREE-DIMENSIONAL SUPERSONIC FLOW 
WITH APPLICATION TO INCLINED 
BODIES OF REVOLUTION 

By John V. Rakich 


Ames Research Center 


SUMMARY 


Characteristics theory for three-dimensional compressible flow is 
reviewed and the compatibility equations are derived in terms of pressure and 
stream angles as dependent variables. The relative merits of bicharacteristic 
and reference plane methods are discussed. A reference plane method is 
developed and demonstrated for pointed and blunted bodies at angle of attack. 

The fundamental complications arising in a three-dimensional method of 
characteristics are that: (1) the compatibility equations contain "cross- 

derivatives M in a noncharacteristic direction; and (2) there is an increased 
need for interpolation to prevent the computed data surfaces from becoming 
skewed. These problems are minimized by using a reference plane method with a 
prescribed and uniform spacing of mesh points. The finite difference mesh 
employed consists of an equal number of points between the body and shock sur- 
faces in each reference plane. Numerical procedures for differentiation and 
interpolation ensure second-order accuracy in terms of mesh spacing. Fourier 
analysis is employed in the circumferential direction and is found to be effec- 
tive in reducing the number of reference planes and computing times. A typi- 
cal mesh consists of 7 planes with 15 points in each plane. The unit 
computation time is about 0.0013 minute per point on an IBM 7094 computer. 

Results for the flow around inclined, circular cones, both blunted and 
pointed, are presented to demonstrate the method described. Predictions of 
surface and pitot pressures are in good agreement with available experimental 
results for a 15° sphere-cone at 10° angle of attack. Bluntness effects on 
shock-layer properties far from the nose are reasonably well predicted. 


INTRODUCTION 


The method of characteristics has been well known for many years and 
excellent theoretical developments can be found in numerous texts and reports 
of which references 1 through 3 give the most complete reviews of multidimen- 
sional theory. However, until recent years there have been relatively few 
practical applications of the methods to three-dimensional flows. This is 
probably due partly to the large number of operations involved in finite dif- 
ference calculations in three dimensions, and partly to the extra degree of 
freedom that results from the existence of characteristic surfaces rather than 



characteristic lines. High-speed computers have made it technically feasible 
to perform calculations for three-dimensional steady and unsteady flows and 
several groups have reported such work. For example, many recent variations 
on the method of characteristics are reported in references 4 through 10, and 
a successful use of a noncharacteristic finite difference method is described 
in reference 11. A number of the proposed characteristics schemes have been 
described and discussed critically in reference 12. 

All the proposed difference schemes for three-dimensional characteristics 
can be placed in one of two broad categories: (1) reference plane methods 

(called semi-characteristic methods in ref. 12) ; and (2) bicharacteristic 
methods. In reference plane methods the characteristic lines are obtained 
from the projection of Mach cones and streamlines onto a prescribed plane. In 
bicharacteristic methods the characteristic lines are the generators of the 
Mach cone and the actual streamlines. Since, in three dimensions, the equa- 
tions written along these characteristic lines are partial differential equa- 
tions, a numerical evaluation of ’’cross-derivatives" off the characteristic 
lines is necessary. (In this sense, three-dimensional characteristic methods 
are similar to noncharacteristic methods for two-dimensional flow.) The prob- 
lem encountered with most bicharacteristic methods arises from the need to 
employ a two-dimensional array of data to evaluate the cross -derivatives and 
to perform the required interpolation in the initial data surface. Further- 
more, if a conventional characteristics mesh is used, the distribution of mesh 
points in the data surfaces tends to become very nonuniform, which causes 
additional complications. For these reasons one is led to reference plane 
methods, where a better control can be maintained on the finite difference 
mesh and where curve fits can be made with respect to a single, variable. 
Reference plane methods have been criticized on the theoretical grounds that 
the initial data may be outside the domain of dependence of the calculated 
point. However, problems related to this criticism have not materialized. In 
fact, the use of such data is precisely what is required by the well-known 
Courant-Friedrichs-Lewey stability condition. 

In the present work, the compatibility equations of characteristic theory 
are derived for both the bicharacteristic and reference plane methods. Prac- 
tical difficulties with the bicharacteristic method are discussed and a finite 
difference scheme based on the reference plane method is proposed. The 
proposed method abandons the usual characteristic mesh and employs uniformly 
spaced points along lines lying in equally spaced meridional planes. Numeri- 
cal interpolation and differentiation are accomplished by means of second- 
degree polynomials in the radial direction and by Fourier analysis in the 
circumferential direction. 

Preliminary results by the present method were described in reference 13 
and compared with calculations by the method of reference 4. Extensive com- 
parisons with experiment were shown in reference 14, establishing the reli- 
ability of the numerical methods. In the present report the flow equations 
and numerical techniques are described in greater detail than was possible in 
references 13 and 14. 
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PRINCIPAL SYMBOLS 


C 1 > c 2^ c 3 
C!*,C 2 * 


x > r > 1 


n,l,t 


u,v,w 


x,r,$ 


speed of sound 
bicharacteristic directions 

characteristic directions in meridional planes 
pressure coefficient 
pitot-pressure coefficient 

unit vectors along x, r, $ 

enthalpy 

total enthalpy 

smoothing constant (eq. (84)) 

coefficient of numerical diffusion term (eq. (84)) 
Mach number 

shock normal and tangent vectors 
pressure 

body nose radius 

projection of streamlines on meridional planes 

streamline coordinates 

velocity components along ©$ 

total velocity V = Vs = Vxj 

axial distance from blunt nose (body axis) 

unit vectors along s, n, t 

cylindrical coordinates 

unit vectors along characteristic coordinates 
unit vectors along 5, n, ? 
angle of attack, deg 


3 


e 

Y 

6 


£ 


ij 


5,n,c 


e 

x 

y 

y* 

v. 

i 

p 

o 

4 


A,B 

b 

j > T 

k 


n 

s 

oo 


transformation matrix (eq. (16)) 
n/M 2 - 1 

specific heat ratio 

shock angle in cross plane (eq. (A6)) 
transformation matrix (eq. (A18)) 
nonorthogonal shock- layer coordinates (fig. 3) 

flow angle from x axis in meridional plane, tan -1 (v/u) (see fig. 1) 

angle between e^ and n 

Mach angle, sin -1 (1/M) 

projection of y on meridional plane 

transformation matrix (eq. (A13)) 

density 

shock angle in meridional plane (eq. (A5)) 
crossflow angle, sin' 1 w/V (see fig. 1) 

azimuthal angle, cylindrical coordinates (see fig. 1), deg 

Subscripts 

initial and new data lines 
body 

indices for radial position of mesh points 
(j = 1, 2, . . . , J - 1, J) (see fig. 5) 

index for circumferential position of mesh point 
(k = 1, 2, . . . , L - 1, L) 

body nose 

shock 

free-stream condition 
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Superscripts 


i index for axial position of mesh point 

v index denoting variables p, 6, <f>, p for v = 1, 2, 3, 4 
* quantity or component referenced to meridional plane 

" unit vector 

-> general vector 

Note: See equation (14) for definition of index notation, 

THEORY 


Basic Equations of Inviscid Flow 

The derivation of the characteristic relations will be carried out in 
this section starting with the equations for inviscid equilibrium flow written 
in their intrinsic form with pressure and flow angles as dependent variables. 
This choice of variables eliminates entropy derivatives from the flow equa- 
tions, thereby simplifying the analysis. Following the development of refer- 
ence 15, the combined continuity and momentum equations are written in vector 
form. 



pV 2 


grad p + div s = 0 


CD 


1 

pV 2 


n 


grad p + t • curl s = 0 


( 2 ) 


— — t • grad p - n • curl s - 0 
pV 2 


(3) 


where s, n, t are orthogonal unit vectors with s tangent to the stream- 
lines. For rotational, inviscid, nonheat-conducting, and nonreacting flow, 
entropy is conserved along streamlines and this requires 

s « grad p - a 2 s * grad p = 0 (4) 
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Conservation of total enthalpy everywhere gives 


h + V 2 /2 = h t = constant 


( 5 ) 


Finally, the equation of state may be written as 


h = h(p,p) (6a) 

or 

a = a(p,p) (6b) 

In the equations above, if n is chosen to lie in meridional planes through 
the body axis, the unit vectors can be written as follows in terms of flow 
deflection angles 0 and <j> (see fig. 1) . 

A 

I* 





(a) Rotation about e^. (b) Rotation about n. 


Figure 1.- Cylindrical coordinates. 


6 


( 7 ) 


s = cos <j> cos 0 e + cos d> sin 0 e + sin <b 

r X Y r Y $ 

n = -sin 0 e + cos 0 e (8) 

x r 

t = -sin d> cos 0 e - sin d) sin 0 e + cos <k (9) 

x r 3> v ' 

Using equations (7) through (9) with the standard vector formulas in cylindri- 
cal coordinates, one can obtain the following from equations (1) through (4): 


B 2 9£ 

P V 2 85 


, 30 

+ COS <b -r— 

y 8n 




8 cj) 

at 


cos <p sin 9 
r 


( 10 ) 


1 

pV 2 


lE. 

8n 


+ 


cos 


30 _ sin 2 <j> cos 0 
3s r 


( 11 ) 


3<j> _ sin <j> sin 0 
py2 r 


( 12 ) 


3p _ _1_ 3p_ 
3s ~ a 2 3s 


(13) 


Characteristics Theory 


Notation . - Characteristics theory is most conveniently developed, 
especially for the multidimensional case, if one uses index notation. There- 
fore, the notation and development of Courant and Friedrichs (ref. 1, p. 75) 
will be followed in this review. Summation convention is used wherein a sum- 
mation symbol is implied by a repeated index. Thus equations (10) to (12) may 
be written simply as 


i 3u V 
a ~ 
yv dx i 


f 

y 


(14) 
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where the indices refer to: 


v dependent variable (p, 6, <j> for v = 1, 2, 3) 

i independent variable (s , n, t for i = 1, 2, 3) 

y equation (eqs . (10) -(12) for y = 1, 2, 3) 

Coordinate transformation.- Equation (14) can be expressed in terms of 
new coordinate directions yj obtained by a rotation of the original coordi- 
nates . Thus one can write 


Xi - ot-p j y j (15) 

where and yj are unit vectors along x^ and yj , respectively, and the 

elements of the transformation matrix are the direction cosines defined as fol- 
lows in terms of the scalar products of x± and y j : 


a 


ij 



X 1 

• yi 

*i 

• y 2 


• y3 

X 2 • 

• yi 

x 2 

• y 2 

x 2 

■ y3 

x 3 

■ yi 

x 3 ' 

■ y 2 

CO 

< X 

■ y 3 


(16) 


With this transformation, equation (14) becomes 


where 



du 

dy i 


f 

y 


b 


3 

yv 


i 

a a. . 
yv ij 


(17) 


( 18 ) 


Char acteristic directions and compatibility equations.- If initial data 
are given for u v on the surface y^ = 0, equation (17) can be solved for 
3u v /3yi in order to generate data on an adjoining surface = Ayi . Plac- 
ing derivatives with respect to y 2 and y 3 on the right side of equation (17) , 
one obtains 


b (1 ^ g 

yv dyi 6 y 


( 19 ) 
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where 


g 


V 


f - b k |Hi 

y yv 3y k 


(k = 2, 3) 


Equation (19) can be considered as a set of algebraic equations for 
3u^/9y^, with a solution by Cramer's rule giving 


3u 

3yi 


yv 


( 20 ) 


where 


yv 


is the determinant of 


b*- 1 ^ and D v 


yv 


is the determinant 


obtained by replacing the vth column of b v J with g 


0 ) 


If 


yv 


yv ~y 

vanishes, then y \ is normal to a characteristic surface and 

the flow equations give no information about derivatives in this direction; 
that is, du v /dyi may be discontinuous. Therefore, the equation of the char- 
acteristic surface is given by 


yv 


( 21 ) 


On the other hand, if equation (21) is satisfied, then the numerator of 
equation (20) must also vanish in order for a solution to exist. Therefore, 


D V = 0 


(v = 1, 2, 3) 


( 22 ) 


gives the so-called compatibility equations. It may be noted from equation 
(19) that the compatibility relations contain one less dimension than the 
original differential equations. For the three-dimensional problem the com- 
patibility equations involve derivatives in two directions. 


Bicharacteristic Method 

The characteristic relations reviewed in the previous section can now be 
specialized to the problem of equilibrium three-dimensional gas flow. Con- 
sider equations (10) through (13) with the supplementary conditions (5) and 
(6). Note that equation (13) is already in characteristic form since it 
involves derivatives in one direction; the streamline, s, is a line across 
which the density gradient, 3p/8n, may be discontinuous. Therefore, one need 
only consider equations (10), (11), and (12) for which the i coefficient 
matrices may be written 
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<$2i cos $ 




B 2 



1 

2 

PV 

1 



«1 


i COS <p 


0 



0 



where 6, . is the Kronecker delta 
ki 

{ 0 k t i 

1 k = i 

The transformed coefficient matrix given by (18) becomes 


b 


(j) 

yv 



B 2 

2 

pV 


1 



1 

2 

pV 


0 l 2 j cos <j> 
aij cos (p 


0 



( 23 ) 


(24) 


Characteristic cone .- The characteristic determinant, equation (21), can 
now be evaluated to give the equation of the characteristic surface in terms 
of the coefficients of the governing differential equations. Using equation 
(24) one obtains 


b^ 

yv 



(25) 


where g 2 = M 2 - 1. Now, from equation (16) the terms in equation (25) are 
recognized as 


“11 = *1 • yi 


«21 = x 2 • yi 


a 3l = £3 • yi 
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which are the components of y\ along 
the x. coordinates. Thus, equation 
(25) describes a cone around the 
or s axis, as indicated in figure 
2(a), making the angle 90 - y with 
the s direction. This cone is nor- 
mal to the characteristic cone. The 
vector yi lies along a generator of 
the normal cone, and the vanishing 
determinant ( 21 ) means that the deriv- 
atives with respect to y\ may be 
discontinuous; that is, the differen- 
tial equations ( 1 ) , ( 2 ) , and (3) 
cannot give any information about 
derivatives in this direction. 

It follows, then, if the y. coordinates are orthogonal, that y 2 and 
y 3 are tangent to the characteristic cone. If y 2 is chosen to lie along a 
generatrix of the cone, then y 2 is called a bicharacteristic direction. 

Compatibili ty equations .- Any one of the three conditions given by 
equation (22) can be used to obtain the compatibility equations. Consider the 
case for v = 1 . 


gl 

ot 2 i cos 4> 

«31 



g 2 

a ll cos ♦ 

0 

= 0 

(26) 

g 3 

0 

an 




Expanding (26) one obtains 

gl " ( a 2l/ a ll)g2 _ ( a 3l/ a ll)g3 = 0 C 27 ) 

where 

81 ’ *1 - “4- (“12 I 77 + “13 § 7 ) - *(“22 §7 * “23 ffj) 

' (“32 “33 |Jj) 

82 * iz - ^2 (“22 |fj * «23 - COS *(«12 * 013 fjj) 

pV \ ' x 


(28) 


(29) 


Bicharacteristic 



Characteristic 

cone 


V 

Normal 

cone 


(a) Characteristic cone and bicharacteristics. 


Figure 2.- Characteristic coordinates. 
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(30) 


- f 3 - -T- (< 

pV v 


“32 


3P 

dy 2 


“3 3 



i$_ + 
3y 2 


“13 


34 > \ 
9 y3/ 


The bicharacteristic direction y 2 can be arbitrarily chosen to lie along any 
ray of the characteristic cone. This means that an infinite number of equa- 
tions can be obtained from (27). However, three equations are sufficient 1 to 
determine the solution for the three dependent variables p, 0 , <£. 


The bicharacteristics can be chosen so as to simplify the compatibility 
equations. First it is noted from equation (16), identifying (x^, x 2 , X 3 ) 
with (s , n, t) , that 




Thus a number of the elements of a. . 
will be zero if y 2 is in the s-n 1 - 1 
plane, and y 2 lies along the t 
axis. In this case (see fig. 2(b)) 
there are two possibilities given by 


sm y 


a . . =1 +cos y 
il 1 


cos y 
±sin y 
0 


(31) 


(b) Bicharacteristics in the s-n plane. 



(c) Bicharacteristics in the s-t plane. 


where the upper sign refers to the 
left-running characteristic C l5 and 
the lower sign to the right-running 
characteristic C 2 . 


Letting 

Yl 

lie in the 

s-t 

plane so that 
axis gives 

y 3 

lies along 

the 

/sin y 


cos y 

°\ 

a. . = | 0 

l 


0 

A 


Figure 2 .- Concluded. 


-COS y 


sin y 


( 32 ) 


where the signs correspond to the bicharacteristic direction C 3 which 
increases with increasing s and t as shown in figure 2(c). 


Redundant schemes which make use of more than three equations are 
discussed in references 8 , 10 , 12 , and 16. 
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Substitution of ou ^ from equations (31) and (32) into equation (27) 

results in three equations along the bicharacteristics Cx, C 2 , and C 3 . 
Using equation (31) for Cx and C 2 , one obtains 


gl ± cot y g 2 = 0 

where y 2 = Cx for the upper sign and y 2 = C 2 for the lower sign. Simi- 
larly, one obtains from equation (32) the following equation for the direction 
C3: 


gx + cot y g 3 = 0 


Straightforward substitution of the g. defined in equations (28) to (30) 
with appropriate elements of ou . from equation (31) or equation (32) yields, 

with some rearrangement, the following compatibility equations: 


6 3P 


PV 


2 3C X 


+ cos 


= ( f i + 6£ 2 - If) si 


m y 


(33) 



3P 

3C 2 


cos 


90 

3C 2 



B£ 2 “ If) Sin y 


(34) 


pV 2 3C 3 3C 3 



cos 


ae\ 

dn) 


sm y 


(35) 


where, from the right side of equations (10) through (12), 

~ cos <f> sin 0 ~ sin 2 <f) cos 0 , ~ sin cj> sin 0 

For two-dimensional flow (<j> = 0) , equations (33) and (34) reduce to the usual 
compatibility equations and equation (35) becomes the streamwise momentum 
equation . 

Fundament al complications . - In comparison with axially symmetric flows, 
equations (33) through (35) have two complicating features which were men- 
tioned earlier. These are (1) the presence of ,f cross -derivatives” 3<j>/3t and 
30/3n on the right side, and (2) the need to perform a two-parameter inter- 
polation for data in the initial data surface. The second problem arises 
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because bicharacteristics through a general mesh point (in a prescribed sur- 
face) do not, in general, pass through mesh points in the initial data 
surface . 2 

Many schemes have been proposed using equations of this form along 
bicharacteristics (see, e.g., refs. 5-8 and 12). However, the programming of 
such methods tends to be cumbersome, and the accuracy of evaluating the cross- 
derivatives is in most cases less than that of the basic calculations. (An 
exception is the recent work reported in reference 10.) These problems are 
minimized if characteristics lying in prescribed reference planes are 
employed. Then the interpolation for initial data and evaluation of cross- 
derivatives can be reduced to a set of one-parameter curve fits with second- 
order accuracy. In the next section, the compatibility equations will be 
derived for characteristic lines which are the projections of bicharacteris- 
tics Ci, C 2> C 3 onto reference planes. 


Reference Plane Method 

In the following development, the flow equations are written in terms of 
two components lying in predetermined reference planes and a third component 
directed out of these planes. For most problems encountered in external aero- 
dynamics it is convenient to specify the reference planes as the axial planes, 
$ = constant, of a cylindrical coordinate system (see fig. 1). The solution 
of problems with axial symmetry is determined by calculation along a single 
plane, but three-dimensional problems require calculation along several 
planes simultaneously. Characteristic theory is employed to calculate the 
flow along each plane. To achieve this, the compatibility equations along 
the projections of the bicharacteristics on the reference planes must be 
derived. The procedure will be to find the pertinent characteristics and 
compatibility relations from the intrinsic momentum equations applicable to 
the reference planes. First, however, it will be necessary to discuss the 
coordinate mesh which will be used to describe the shock layer. 

Finite difference mes h.- The cylindrical coordinate system used to 
expand the vector relations in equations (1) through (3), and to define the 
reference planes, is not ideal from the computational standpoint. This stems 
from the fact that special treatment would be required for bodies with non- 
circular cross sections and, more importantly, for shock surfaces. One is 
therefore led to a finite difference mesh which divides the shock layer into 
a number of annular rings which include both the shock and body surfaces, as 
shown in figure 3. The resulting mesh points are connected by a nonorthogonal 
£, n, C system of coordinates: £ and n lie in the reference plane with n 
usually normal to the body surface; £ is directed out of but not generally 
normal to the reference plane. 


2 Conversely, bicharacteristics through known points on a plane initial 
data surface will, in general, intersect at points lying in a nonplanar sur- 
face. The subsequent data surfaces become increasingly distorted as a compu- 
tation proceeds away from the initial data plane. 
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Figure 3.- Nonorthogonal shock- layer coordinates. 


A related system is one in which 
5 is replaced locally by the projec- 
tion, s* , of streamlines onto the 
reference planes. It is this s*, n, 
C system which is used to express 
the intrinsic flow equations (10) 
through (13) in a form needed for the 
present reference plane analysis . 

The unit vectors, 5L = (s , n, t) , in 

streamline coordinates are related to 
the new system by direction cosines, 
e * . , defined by 


x. = e * . z .* (36) 

l 13 3 


where z_.* = (s*, r \ , £) . 

Written in terms of its components, equation (36) gives 


s = 


e* x s* + e* 0 n + e* 0 c 


'12" T c i: 


(37) 


n = 


e* x s* + e* 2 n + e* 3 c 


(38) 


t = 


e 31 s* + e* 9 n + e* q £ 


"32" ^ 33 " 


(39) 


Appendix A describes how is calculated in terms of the body and 

shock-wave shapes. The detailed expressions for the direction cosines are 
not needed for the development of this section but it should be noted that, 
since the shock shape itself is obtained from a solution of the problem (i.e., 
a direct as opposed to an inverse approach) , is not known beforehand for 

the entire flow. However, in a locally supersonic region, where the shock can 
be calculated step by step, e*. can always be determined as the calculation 

proceeds downstream from an initial data line. 

Reference plane equations .- Recalling the rules for a directional 
derivative, equations (37) and (39) yield the following expressions for 
derivatives in the s and t directions. 
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3 _„* 3 3 3 

as - = e n + e i2 77 + £l3 7c (40) 

3 _ * 3 *3 *3 

dt ~ e 31 + £ 32 37 + £ 33 77 (41) 


These differentiation rules allow one to write the intrinsic flow 
equations in terms of the desired planar components. Substituting equations 
(40) and (41) into equations (10) through (13) and regrouping terms, one 
obtains 

* e 2 a P . ___ , ae 


£ n — 7“^ + cos 


pv as 


= -F * 

an x l 


(42) 


pV 


i 9p * , ae 

~ 77 + £ n cos $ 


= f. 


as 


(43) 


as 


= f 3 


(44) 


9p _ f * 
“ ft * 

3 s 


(45) 


The left sides of these equations contain derivatives in the reference planes 
and the remaining terms are all lumped into the f^*, which are given by 


fi* - - 


: 2 = 


COS <f) 

sin 

e 
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* 9 p\ A 
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CM 
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31 * 12 3f| ^ 3 9 £ 
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(46) 
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The f . * thus expressed contain derivatives along the n and £ coordinate 

directions and, as a result, can be easily evaluated with standard one- 
parameter differentiation routines. It will also be seen, when the finite 
difference scheme is described, that the derivatives 3cj)/3s* appearing in 
f^* and 3p/3s* appearing in f 3 * and f 4 * are easily treated. 

It is important to note that, although written in a simplified form, 
equations (42) through (45) contain no additional approximations beyond those 
in the original equations. Only the usual assumptions pertaining to viscosity 
and heat conduction are made. The bracketed [ ] terms in the f . * all vanish 
for axially symmetric flows (i.e., for <$> = 0) , and the equations reduce to 
the familiar intrinsic equations for zero angle of attack. 

Characteristic direc tions . - Characteristics theory, as outlined in a 
previous section, can now be applied to equations (42) through (45). It is 
first observed that equations (44) and (45) are already in the desired char- 
acteristic form. These equations give no information about the normal deriv- 
ative 3/3n; therefore, s* is a characteristic direction for c|> and p. This 
is in contrast to equations (42) and (43), which can be combined to give 
derivatives in different directions. Writing the latter in the form of 
equation (14), one has 


= f * 
3x^* M 




(50) 


Expressed in terms of new coordinates, y.*, equation (50) becomes (cf. eqs. 
(17) and (24)): 1 


where 



(51) 


and where 


a. . 


are the elements of the transformation matrix relating x. 


and y . The characteristic directions for equation (51) are obtained from 
the determinant 
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b * o ) 

pv 


= 0 


which yields 



( 52 ) 


Writing 
for this 


a . . 
case 


in terms of a rotation angle as in equation (31) , one obtains 


a . . 
ij 


/ sin y* 

V +cos y* 




where y* is the angle between the 
streamline and the characteristics 
in the reference plane as shown in 
figure 4. Equation (52) shows that 
y* is related to the Mach angle 
according to 


Figure 4.- Characteristic directions for the 
reference planes. 


cot p* = e* x cot p (53) 


Compatibility equation s.- The compatibility equations applicable to the 
directions Cj* and C 2 * can be obtained in the way previously described for 
the bicharacteristics. For the present case, one has in place of equation (27) 


where 


and 


“21 

81* - SIT 82* = 0 


81* - fi* - hip *2™. 


9y 2 


9y 2 


CO 


82* * 82 * - s™ i- 


,C2) 


9y 2 


(54) 


Algebraic details are straightforward and need not be repeated. The resulting 
compatibility equations are: 
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+ COS (j> 


90 


= (fi* + gf 2 *)sin y* 


(55) 


B 3P 
2 * 
pV 3C i 


9C i 


* 


-^2 3p - COS <j> = (fi* - Bf 2 *)sin y* (56) 

pV 3C 2 * 3C 2 * 

Equations (55) and (56) , together with equations (44) and (45) , are the final 
set of differential equations programmed for numerical calculations . They are 
supplemented by the energy equation (eq. (5)) 

h + V 2 /2 = H 

and the equations of state ((6a) and (6b)) 

h = h(p,p) 


and 


a = a(p,p) 

Details of the numerical methods used in the computer program are described 
in the next section. 


NUMERICAL METHODS 


The general theory of characteristics was developed in the previous 
section, where it was shown that the three-dimensional problem can be reduced 
to an equivalent two-dimensional form. A numerical solution of the problem 
can then be accomplished in the usual way by numerically evaluating deriva- 
tives in one direction in order to calculate a step forward in the second 
direction. The problem is analogous to the numerical solution of two- 
dimensional hyperbolic equations by noncharacteristics methods. 

Therefore, in formulating a practical method for calculating three- 
dimensional flow, the numerical differentiation process is of primary impor- 
tance. If the differentiation is to be performed efficiently and accurately 
(i.e., at least second order in a typical mesh dimension), the mesh points 
should be constrained to lie along simple coordinate lines. Secondly, the 
boundary calculations are simplified if the coordinates lie on the shock and 
body surfaces. These considerations suggest the shock- layer-oriented, non- 
orthogonal coordinates shown in figure 3. The resulting computational proce- 
dure is simplified significantly compared with previously proposed three- 
dimehsional characteristics methods (see, e.g., ref. 12). 
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Presented next are the difference equations, computational logic, and 
numerical differentiation procedures developed on the basis of this nonorthog- 
onal (5, n, c) mesh. The problem may be stated as follows: 

Given data on an initial n - C surface at 5 0 where the flow is 
locally supersonic, it is required to generate, by means of the 
flow equations and boundary conditions, new data on the adjacent 
n - £ surface at 5 0 + 


Finite Difference Equations 


Figure 5 shows the mesh point arrangement for a typical reference plane 



and (i) the initial data line at 5^ 
and new data line at £g. Let the 
subscripts j and k denote the radial 
and circumferential positions of the 
mesh points. However, the subscript 
k will be omitted for brevity in 
writing the difference equations. 

The method adopted 3 to calculate 
conditions at a typical mesh point 
i,j on makes use of interpolated 

data at the points of intersection on 
5^ of the characteristics through 
point i,j. A three-point Lagrange 
interpolation is employed to determine 
the necessary data from known condi- 
tions at neighboring mesh points . 


Three intersections are required for each field point on To identify 

these points the convention is adopted whereby the subscript t represents 
the intersection with 5 A of the streamline projection s* through point 
i,j, and the subscripts t-1 and t+1 represent intersections of characteris- 
tics Ci* and C 2 *, respectively (see fig. 5). 

With this convention the compatibi lit) 
written in the following finite difference 

- pi:!)* ».(•/ - 


M p / -Pm)- B *( 9 3 l - 


3 This method is called the Hartree method (ref. 17) and also the inverse 
method (ref, 12); it was used by Katskova and Chushkin (ref. 9). 


r equations (55) and (56) can be 
form: 


e^j) = ? 1 ACr 


(57) 


) 1 1 ) = 
T + 1 / 


(58) 


e :. = p 2 AC 2 * 
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Equations (44) and (45) are similarly written as 


'J’j 1 - = F 3 As* ( 59 ) 

and 

Pj 1 - P^” 1 = F 4 As* (60) 


The system of equations is completed by the energy and state equations, (5) 
and (6). They apply to all of the field points - that is, for the index j 
running from 2 up to J-l. For the body point, j-l, equation (57) is replaced 
by the equation of the body. 


rj 1 = f (x ! 1 , § k ) 


which permits the calculation of e^ 1 by means of equations derived in 
appendix C. It is shown there that 


where 



(61a) 


For bodies of revolution equation (61a) reduces to 


e l 


i 


tan 1 


a *! 1 

9x 


(61b) 


m 

At the shock, j = J, equations (58), (59), and (60) are replaced by the 
oblique shock equations. The jump conditions for a general three-dimensional 
shock surface are developed in appendix B and can be written in the following 
functional form for fixed free-stream conditions: 





(62) 
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where u V represents p, 0, <j), p for v = 1, 2, 3, 4. Here a and 6 are 
the shock-wave angles in the planes $ = constant and x = constant, respec- 
tively, and a is the angle of attack. Equation (62) depends directly on the 
unknown shock angle a 1 , and indirectly on the parameters 6k, and a. The 
shock angle 6^ is determined by numerical differentiation of shock coordi- 
nates as described in appendix A. The four equations obtained from (62) for 
v = 1 through 4, together with equation (57), are sufficient to determine the 
shock angle a 1 . 


Computational Procedure 

Local iteration.- The difference equations are solved by a standard Euler 
predictor -^corrector method in which the coefficients are treated as constants 
locally and equal to their average value over the step. An initial guess, say 

(u v )^ = (u V )j 1 , is used to start an iterative procedure by which corrected 

values of (u v )^ are calculated using coefficients evaluated with the previous 

predictions. The iteration is continued until the pressure repeats to a 
specified accuracy. For typical mesh points three correctors reduce the rela- 
tive error to less than lxlO" 5 . It is shown in standard texts that the iter- 
ated result has a truncation error of the order of the step size cubed. 

In this method the coefficients Ay, By, and Fy in equations (57) through 
(60) are averaged along the characteristic direction appropriate to each 
equation. For example, the coefficient Fi in difference equation (57) repre- 
sents the average of the right side of differential equation (55) , taken along 
the characteristic Ci*, and is written as 

F i = l/2[(fi* + Bf 2 *D]j + 1/2 [ (f i* + (63) 

in the present notation. Similarly, averages are evaluated using data at 
point t+ 1 for equation (58) and at point x for equations (59) and (60) . 

Global iteration.- The set of equations (57) through (60) are solved 
successively on L reference planes, \ (k = 1 , 2 , , . . , L) . The differ- 
ence equations governing the flow along various reference planes are coupled 
by cross-derivatives which are included in the functions Fy appearing on the 
right side of each equation and by the shock angle 6k appearing in equa- 
tion (62) . In order to solve these equations in an explicit manner, it is 
therefore necessary first to approximate Fy and 6^ with derivatives evalu- 
ated on 5 a- Then, using calculated values on all of the L planes to 
evaluate cross-derivatives on one obtains the next approximation to Fy 

and 6 k and the entire process can be repeated. This is referred to as a 
global iteration, in contrast to the point-by-point iteration employed in the 
local solution of the difference equations. 
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Since computing times are generally large in three-dimensional problems 
(see next section), test cases were run to determine the need for global 
iteration. Table I shows that excellent results are obtained without itera- 
tion, that is, using derivatives evaluated on It is therefore expected 

that this iteration would not be required in most problems. 


Step Size 

Since the mesh points are not constrained to follow characteristic lines 
with the present method, the size of a forward step, AC = can be 

arbitrary to some extent. The only limitation is that, in order to have a 
stable numerical process, the step should not exceed a certain maximum deter- 
mined by the region of influence of the initial data. This stability condi- 
tion will make AC depend on the radial step An, smaller radial steps 
requiring smaller forward steps. The lateral step Ac does not affect the 
stability directly, although it does, of course, affect the accuracy. The 
effect of step size on accuracy is discussed after the following paragraphs 
on stability condition. 


Stability co ndition .- Although the 
analysis of numerical stability has not 
been performed for the general nonlin- 
ear equations, the stability criterion 
for linear hyperbolic equations (see, 
e.g., ref. 17) is apparently sufficient 
for the nonlinear equations . This is 
the well-known C-F-L (Courant, 
Friedrichs, Lewey) condition which 
essentially states that the domain of 
dependence of the calculated point must 
be included in the initial data. This 
means that the characteristic C}* 
through point (i,2) in figure 6 must 
pass through or above point (i-1,1). 
Similarly, the characteristic C 2 * 
through point (i,J-l) must fall through 
or below point (i-l,J). 

Strictly followed, the C-F-L condition would require testing every 
shock and body point to determine the maximum allowable step size, AC m , which 
would insure stability. However, experience indicates that the condition is 
not overly restrictive in the sense that steps slightly larger than AC m do 
not usually cause a violent instability. Therefore, it is adequate to test 
only at the body and, for bodies with circular cross section, on the windward 
side where AC m is likely to be smallest, due to the low Mach number. The 
step size is then taken slightly less than the maximum; a value of 
AC = 0.8 AC m works well in most problems. 

Accuracy and computi ng time . - The accuracy of a numerical computation is 
usually estimated by comparing results obtained with various mesh spacings . 



Figure 6.- Step size limitation. 
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Since the exact solution is usually unknown, one can only compare with the 
results of the finest mesh and observe how the error decreases. If the trun- 
cation errors are second order in terms of the mesh size, then halving the 
mesh size should reduce the error by one-quarter. 

To test the accuracy of the present method, calculations were performed 
from x/R n = 2 to x/R n = 3, using 3, 5, and 7 planes and 5, 10, and 15 points 
along each plane. In each case the same starting conditions were used at 
x/R n = 2. The results of this study are shown in table II. Table 11(a) pre- 
sents the shock angles and the surface pressure on the $ = 90° plane at 
x/R n = 3. These results show that the method is of second-order accuracy. 
Also, it should be noted that the results with k = 5 and 7 agree very well 
for J fixed (moving horizontally in the table). The good accuracy with only 
a few planes is attributed to the use of trigonometric analysis for the 
cross -derivatives . 

The computing time naturally increases as the mesh is refined; this is 
illustrated in table 11(b), which lists the total number of points computed, 
the total execute time, and the time per point in minutes. The calculation 
was performed on an IBM 7094 Model l^in FORTRAN IV (version 13 IBJOB proces- 
sor) . A unit time of about 0.13x10 min per point was obtained with the 
finest mesh. The unit time increases as the number of points decreases, prob- 
ably because of fixed input/output times. The actual and unit times are 
almost doubled when one global iteration is made at each station. 


Numerical Differentiation 

Discussed next is the problem of evaluating cross-derivatives appearing in 
the equations. Partial derivatives in three directions appear in the functions 
f±* defined by equations (46) through (49). They are of the form 8/9s*, 
d/dr ] , and 9/3?, in the stream, radial, and circumferential directions. Special 
treatment is given to the circumferential derivative, following the discussion 
of the first two. 

Stream wise and radial derivatives.- These derivatives are taken together 
since they are both evaluated by the standard polynomial approach. The main 
idea is to employ a differentiation formula consistent with the accuracy of the 
basic calculation. 

For the streamwise derivative the approximation 

f v. i , v. i- 1 

, tU - (U K 

9s* As* 


(64) 


is clearly equivalent to the form of the difference equations employed. 
Equation (64) represents the derivative of u v at the midpoint of the inter- 

2 

val, As*/2, with an error of the order As* . This derivative is evaluated at 
each step in the local iteration process previously described. 


24 



Because of the equally spaced coordinate mesh presently employed, the 
radial derivative is similarly determined to the same accuracy by the central 
difference 


, v>i , v.i 

.v (u ) . - (u ) . 

3u J + 1 3-1 


9n 


2 An 


(65) 


Equation (65) approximates the derivative at point (i,j) with an error of 
order An 2 . The average radial derivative at the midpoint (i+l/2,j) between 
and can be obtained by means of the global iteration procedure 
previously described. 

Standard end-point formulas of equivalent accuracy are used at the body 
and shock where central differences are not possible. These need not be 
written out, as they can be found in many books (e.g., ref. 18). 

Circumferential derivatives . - In the aerodynamics of bodies it is well 
known, from linearized and perturbation theories as well as from experiment, 
that the solution of most problems can be represented by a trigonometric 
series in the meridional angle $. When information such as this is avail- 
able it should be possible, by choice of an appropriate functional form, to 
improve a numerical process. For example, with data known to have a nearly 
cosine variation, it is clear that fewer sample points are required to 
approximate the data with a cosine series than with a polynomial. A 
Fourier-series approximation is therefore used to evaluate derivatives with 
respect to £. Symmetry conditions, which usually arise at $ = 0 and $ = tt, 
are easily satisfied in this way. This technique makes it possible to cal- 
culate with fewer reference planes, thereby increasing the computational 
efficiency (see ref. 13) . 

For the present application it is necessary to determine Fourier 
approximations for pressure p, flow angle 0, crossflow angle c|>, and 
density p from their values on L planes \ (k = 1, 2, . . . , L) with 
$1=0 and $l = tt. Symmetry conditions for three variables (represented by 
u k v (v = 1, 2, 4)) are satisfied by a cosine series of the form 




V 


L-l 

2^a n V cos n 4> k (v = 1, 2, 4) 

n=o 


( 66 ) 


For the crossflow angle (v 
series 


3), which is zero at $ = 0 and $ = tt, a sine 


L-l 



sin n 

k 


(67) 


is necessary. 
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When equally spaced meridional planes 0^ = 7i(k - 1)/(L - 1), 

(k = 1, 2, . . . , L) are used,, the calculation of the Fourier coefficients 
is particularly simple because the usual orthogonality conditions are 
exactly satisfied by a finite sum (see, e.g., ref. 18 or 19). In this case 
the coefficients are given by 


and 


v 

a 

n 




+ u T 


cos nTT 


) 



(n = 0, 1, . . L - 1) 


( 68 ) 


i 


b n = F sin n \ (n = 1, 2, . . . , L - 1) 

k=2‘ 


(69) 


Derivatives with respect to the angle 0 are obtained by differentiating the 
Fourier series, equations (66) and (67), to give 



dO 


L-l 


-na sin n 0, 
n * 


n=l 


(70) 


and 


, L-l 

v 1 

v , nb n C0S r \ (71) 

n=i 


The desired derivatives in terms of distance along the £ direction are 
obtained from equations (70) and (71) according to 


3 

as 



s 


(72) 

if 


where it is understood that the derivative on the right side of (72) is 
evaluated using data on the £ coordinate. The scale factor g^ relates 
the distance along c corresponding to an incremental change in 0, 
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and can be written in terms of the direction cosines developed in appendix A. 
From these coordinate relations one has 


r(d$/dO = e $ • i = v 33 

and therefore, 

g ? = (r/v 33 ) 


(73) 


4 


* 


4 


4 


Starting Data 

The method of characteristics, being a method for initial value 
problems, requires starting data which are usually obtained from boundary 
conditions or from an initial solution obtained by other methods. For pres- 
ent applications two types of starting conditions are of particular interest. 
These are (1) the spherically blunted body with the sonic line located on the 
spherical nose, and (2) pointed bodies which can be approximated by a cone in 
some small region near the tip. 



Sphere . - Since a sphere does not have a 
preferred orientation, the flow remains symmet- 
ric around the wind axis for any angle of attack 
Therefore, axisymmetric blunt-body solutions 
obtained, for example, by the inverse method of 
reference 20 can be used to provide initial con- 
ditions. It is necessary, however, to generate 
these axisymmetric starting data on an initial 
n-£ surface which is defined in a body axis 
system; the data will not be symmetric with 
respect to the body axis . 




Figure 7.- Starting data for a 
spherical nose. 


The characteristics program is set up to 
generate body axis data from wind axis data 
obtained from a blunt-body solution. Given 
these wind axis data on one body normal (see 
fig. 7(a)) where the flow is supersonic, say 
M > 1.05, the characteristics calculation is 
performed in a wind axis system and for a = 0 
to the position locating the body normals 

nk * • (Primes denote wind axes.) This is illus- 
trated in figure 7(a). Since the flow is axi- 
symmetric in terms of wind axes, the normals 
may be placed at an arbitrary circumferential 
position. The values of x^* are chosen so 
that the normals nk f match the ring of normals 
rik emanating from the sphere at (x 0 ,r 0 ) in 
terms of body axes (fig. 7(b)). They are 
related to the meridional angle and the 

angle of attack a by 

x k 1 ~ **n + ( x o ~ R n) cos a - r Q sin a cos $k 

(74) 

where Rjj is the radius of the spherical nose. 
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Scalar data derived in this way on the body normals 
(k = 1, 2, . . . , L) are ready for use in the general three-dimensional cal- 
culation. However, the flow angles, 6 and <j>, must be recalculated in terms of 
body axes by means of the transformation relationships for velocity components. 
By the definition of 0 and <J> (fig. 1) . 


0 = sin 


v 

V cos <f> 


(75) 


4> = sin 1 y (76) 

The magnitude of the velocity is unchanged by the coordinate rotation, so 
that 


V = V' (77) 

From reference 15, the velocity components v, w can be expressed as fol- 
lows in terms of wind axis components : 


v = (u 1 sin a 

+ v T cos a 

cos $ f )cos 0 + v ! 

sin 

O’ 

' sin 

$ 

(78) 

w = v' sin $' 

cos $ - (u T 

sin a + v T cos a 

cos 


1 ) sin 

$ 

(79) 


where $ ' is obtained from 


cot $ sin $' - cos a cos $' 


<*' - V . 

r— sin a = 0 

r * 


Equations (75) through (80) determine 6 and <j>, relative to body- axis 
coordinates, from u T and v T calculated in the wind-axis system. 


(80) 


Cone.- Conical solutions are defined as those which are independent of 
g - that is, all derivatives with respect to £ vanish. Conical flows can 
be obtained in two ways: (1) by solving the boundary value problem for the 

reduced equations with 3/3£ = 0 (see, e.g., ref. 21), or (2) by the asymp- 
totic solution of three-dimensional initial value problem with a conical body 
(see ref. 11 or 22) . The latter approach is presently taken since it fits 
the general computation scheme with little change. 

Approximate initial data for a cone are specified and the calculation 
downstream is carried on until the conical flow condition is met to a speci- 
fied accuracy.^ This approach has the disadvantage of being generally more 
time-consuming than the boundary value method, but it avoids many difficul- 
ties inherent in boundary value problems. The situation is quite similar to 
the calculation of the supersonic blunt-body problem by an asymptotic 


K 


i 


4 It is clear physically that the conical condition must be obtained 
far downstream from the approximate starting condition. 
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unsteady calculation. Approximate starting conditions are obtained either 
from the perturbation solution for cones at small angle of attack (ref. 23) 
or from a previous three-dimensional solution for a different Mach number or 
cone angle. 


Two modifications to the general computational procedure are required 
to obtain the correct conical flow solutions for circular cones. First, it 
was found necessary to account for the Ferri vortical layer (ref. 24) by 
setting the entropy at all body points, except at the leeward plane of 
symmetry, equal to the entropy at the windward shock point. The entropy at 
the leeward body station, which is a singular point, must equal the entropy 
at the leeward shock point. 5 Thus the following conditions apply to circular 
cones : . 


1 ,k 


= S 


J,L 


(k = 2, 3, 


. , L) 


and 


(81) 


S 

1,1 



Secondly, as a result of the entropy singularity at the leeward body point, 
the circumferential density derivative there is also singular. Therefore, 
data from the singular point must be excluded in the numerical differentiation 
for 9p/3£. 


Smoothing 

Under certain conditions where the density or crossflow angle develop 
large radial gradients, the numerical calculation for these quantities appears 
to be unstable. This situation can arise in the region of the vortical layer 
on pointed cones and in the so-called entropy layer over blunted cones. A 
stabilizing difference scheme similar to that known as the Lax, or Lax- 
Wendroff method (refs. 26 and 27) is therefore used in these cases. 

Equations (59) and (60) are modified for this purpose according to the 
following difference approximation: 


, v, 1 
CU ) 3 


, v.i-l 
Cu ) 


/ 


F 

v 


k Aq 3 
+ 2 ^j 




(82) 


where v = 3 or 4 represents <p or p , respectively. The second term in the 
bracket is proportional to the second derivative of u v , so that the 


5 Additional singular points arise in more general conical flows (see, 
e.g., ref. 24 or ref. 25). 
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difference equation (82) approximates the differential equation 6 



F 

v 


+ K 


aV 

3n 2 


( 83 ) 


where 


y _ k An 3 
- 2^ AC 


(84) 


and r) is the value of n at the shock Hj = (J - l)An. 

The additional term in equation (83) represents a diffusion process and 
is the stabilizing term in the type of difference scheme given by equa- 
tion (82) . The arbitrary constant k in equation (84) is included to allow 
control of the diffusion term. If k = 0, equation (83) reduces to equa- 
tion (S9) or equation (60). For k = nj/An, equation (83) is essentially the 
same as the Lax scheme (ref. 26) . 

It is known (this is demonstrated below) that the Lax difference scheme 
provides too much diffusion; therefore, it is desirable to choose k between 
0 and nj/An in order to provide numerical stability without undue loss of 
accuracy. In most applications k has been set equal to 1, and this condi- 
tion is presently termed a second-order Lax smoothing because in this case 

K = 2n A \ g An 2 = 0(An 2 ) (85a) 

J 


Thus, the second-order difference scheme approaches the differential equation 
like An 2 as the mesh is refined, whereas the Lax method approaches linearly 
in An. For large Mach numbers, an order- of-magnitude analysis reveals that 



which gives K a value of about 1/1000 for a typical mesh. 


(85b) 


6 This argument is not strictly rigorous, as pointed out in reference 27, 
if the dissipative term is of the same order as the truncation error of the 
other terms. Nevertheless, the analogy is qualitatively correct, especially 
when the step size is not small. Also, when the integration is carried over 
large distances the truncation error might tend to be random, while the 
dissipative term is always additive. 
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Figure 8.- Effect of smoothing on density 
distribution at x/R r = 40; 15° sphere- 
cone, - 10, a = 0. 


The effect of this difference approximation is illustrated in figure 8 
for a 15° sphere-cone at zero angle of attack. The figure shows the density 
distribution in the shock layer at a station 40 nose radii downstream from 
the stagnation point. 

Calculations by the present method are compared with the method of 
reference 28 which accurately determines the density distribution in axisym- 
metric flow by the use of a stream function. Present results with second- 
order smoothing are shown in figure 8 by the broken lines for 15 and 20 mesh 
points. It is seen that these curves differ from the more exact result only 
where the second derivative, 3 2 p/8n 2 , is large, and that the error decreases as 
the mesh is refined. Results obtained without smoothing and with 15 mesh 
points are shown by the symbols. The calculations appear to be unstable 
without smoothing. 

It is emphasized that the need for smoothing arises only when the 
computational mesh is coarse relative to the details of the flow field. Exper- 
ience indicates that the smoothing can be eliminated if the mesh is refined 
sufficiently. This, however, is not practical in three-dimensional problems 
where computation time and computer storage limitations are important factors , 
Therefore the proposed smoothing scheme provides a means for obtaining meaning- 
ful results with a relatively coarse mesh of points. 

Disconti nuous derivatives .- The possibility of discontinuities in 
supersonic flows presents a severe test for a numerical calculation. With the 
present numerical method, characteristic lines are not followed and some smear- 
ing of the discontinuities will occur because of the interpolation for data on 
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(a) Leeward plane, $ = 0 



(b) Windward plane, $ = 180° 


Figure 9.- Pressure distribution on surface nor- 
mals, n; 30° sphere-cone, = 10, a = 5°. 

should have a discontinuous slope there, 
are only slightly rounded by the quadratic interpolation presently employed 


the initial data line. This problem 
was therefore investigated and it was 
found that the proposed method does 
provide an accurate resolution of known 
flow details obtained with a standard 
method of characteristics for axisym- 
metric flow (ref. 29) . 

Discontinuous derivatives arise 
primarily on the body surface at points 
where the surface curvature is not 
analytic. As an example, figure 9 
shows a map of the pressure in the 
vicinity of the sphere-cone juncture 
for a blunted cone at 5° angle of 
attack. Each line represents the pres- 
sure variation along a body normal, n , 
between the body and the shock wave. 

The origin for each line is displaced 
in proportion to the x coordinate of 
the body point so that the circular 
symbols represent the actual surface- 
pressure distribution. The surface 
pressure decreases uniformly on the 
spherical nose until the sphere-cone 
juncture is reached, at x^ = 0.5, 
where the pressure gradient, 3p/3s, 
changes abruptly. The discontinuity 
from that point moves out into the flow 
field along a Mach wave. (The approxi- 
mate position of the wave is indicated 
by the arrows.) The pressure is nearly 
constant behind the wave and varies 
according to the blunt-nose flow ahead 
of the wave. The curves theoretically 
As seen in the figure, the curves 


RESULTS AND DISCUSSION 


Many numerical solutions have been obtained with the described method of 
characteristics, and these results compare favorably with other published 
works (see, ref. 13). Presented in this section are the details of a typical 
calculation for a 15° sphere-cone at 10° angle of attack and a Mach number 
of 10. This is considered to be sufficiently nonlinear to provide a good 
test of the theory. The dominant three-dimensional features of the flow are 
illustrated and compared with predictions of perturbation methods for small 
angles of attack. Bluntness effects are illustrated by comparison with a 
pointed cone solution, which is described first. 
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Pointed Cone 


The solution for a pointed 15° cone at a Mach number of 10.6 and 10° 
angle of attack is presented in table III. Listed for each plane, 

0 = constant, are the shock angles, a and 6. Below the shock angles are 
tabulated the x,r coordinates of each mesh point running from the shock to 
the body and the corresponding flow variables at those points, as labeled. 
The quantity M* is the component of Mach number defined by the Mach angle 
y* in equation (53) . Remaining variables are defined in the list of 
symbols . 


This conical solution was obtained with a coordinate mesh consisting of 
9 meridional planes and with 11 points on each plane. Second-order smooth- 
ing was used on the density p and the crossflow angle cj> (smooth constant 
k = 1) . Initial data for the present case were obtained from a previous 
solution at M = 7 which agreed with the tabulated results of reference 11, 
The free-stream Mach number was changed from 7 to 10 and the computation was 
carried downstream until the flow relaxed to the new conical solution. This 
was accomplished in a number of stages, with each stage consisting of a com- 
putation from x = 0 . 8 to 1.0, using the output of the previous stage as 
initial data. 


Final Starting 

data data 

M ro = 10.6 M ® =7 



Figure 10.- Relaxation of conical solution after 
N stages of calculation from x = 0.8 to 
x = 1.0; 15° cone, M = 10.6, a = 10°. 


The accuracy of the solution is 
indicated by figure 10 which shows 
the shock angle in three meridional 
planes as a function of the reciprocal 
of the number of stages. It is seen 
that the approximate starting condi- 
tions cause the shock angles to change 
abruptly in the first stage with a 
slow decay to a limiting value as 1/N 
tends to zero. Conditions on the lee 
side are slowest to approach a limit. 
The shock angle at $ = 0 repeated to 
an accuracy of Aa/a = 0.45x10” 3 in 
the last stage of calculation and 
repeated to an accuracy of lxl0~ 5 
during the last step of the tenth 
stage . 


The main features of this conical flow will be illustrated in conjunction 
with the blunted cone results described next. 


Spherically Blunted Cone 

The calculation for a spherically blunted cone with a 15° semi vertex 
angle in air (y = 1.4) at M^ = 10 used 7 meridional planes with 15 points in 
each plane. Second-order smoothing, k = 1, was employed for the 10° angle- 
of-attack solution when x/^ > 10; no smoothing was necessary for the other 
solutions presented. 
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, Shock surface 





<£>, deg 


Figure 12.- Shock trace on body-normal surfaces; 
15° sphere-cone, = 10, a = 10°. 


The shock profiles in three meridional planes are shown in figure 11 for 
a cone length of 20 nose radii. Also shown on the figure are the pointed- 
cone shock positions which the blunt-cone shock must approach as x/R n gets 
large. The small difference in the Mach number of the two solutions of 10 
to 10.6 does not significantly affect the comparison. It is observed that 
the shock quickly approaches the pointed-cone shock for $ = 90° and $ = 180°, 
while it does not appear to approach the pointed limit on the lee side, 

$ = 0° . This is further illustrated in figure 12, which shows the circum- 
ferential shock shape for several axial stations. The slow decay on the lee 
side is evidenced by the hump in the profile near $ = 0. Note that the x 



coordinate is not constant on the shock 
traces shown, because the shock radial 
position is measured at its intersection 
with a cone normal to the body surface 
(see sketch in fig. 12). 

The axial distribution of shock 



angle is shown in figure 13. On the 
windward plane, $ = 180°, the angle 
reaches a minimum value of 13.744° at 
x/Rjj =3.89. On the leeward plane the 
angle has not reached a minimum by 
x/R n = 20 and is still well above the 
pointed- cone value at that point. 

In a three-dimensional flow the 
streamlines generally flow across the 
meridional planes as a result of the 
asymmetrical pressure distribution. 


Rn 

Figure 13.- Shock angles in meridional planes; 
15° sphere-cone, M w = 10, a = 10°. 


This crossflow is described in terms of 
the angle <j> which is related to the 
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Figure 14.- Axial distribution of surface cross- 
flow angle; 15° spherc-cone, = 10, 
a * 10°, <J> = 90°. 

familiar crossflow velocity according 
to 4> = sin -1 (w/V) . Figure 14 shows 
the axial distribution of crossflow 
angle along the side meridian 0 = 90°. 
The crossflow angle is a minimum at the 
sphere-cone juncture and then rises to 
a maximum value about 2.4 times the 
angle of attack. Thus, the surface 
upwash for blunted bodies can be greater 
than the maximum value of 2a given by 
s lender-body theory, which applies to 
pointed bodies at low supersonic speeds. 
It is interesting also to compare the 
result of the linearized characteristics 



O 20 40 60 80 100 120 140 160 180 

<£> deg 


Figure 15.- Circumferential distribution of 
surface crossflow angle; 15° sphere-cone, 
Nt * 10, a - 10°. 


method (ref. 15) which is in good agreement near the nose. Agreement extends 
to about x/R n = 10, where the linear method starts to break down. 


The reason for the breakdown is evident in figure 15 which shows the 
circumferential variation of 4> for a = 10° at various axial stations. Near 
the nose the variation is nearly sinusoidal, as assumed in the linearized 
method. However, the exact three-dimensional calculations show a large 
deviation from the sinusoidal variation beyond x/R n = 10. 

Figure 16 shows the variation of crossflow angle normal to the body. 

The curve for X ^/R n = 20 approaches the pointed- cone solution near the 
shock, n/n s - 1, but deviates by a large amount near the body. This is due 
to the low density of the flow in the entropy layer which is generated near 
the body surface by the blunt nose. The flow in this region has less momen- 
tum that that away from the surface, and is therefore turned more by the cir- 
cumferential pressure gradient. The situation is analogous to boundary- layer 
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Figure 16.- Radial distribution of crossflow 
angle; 15° sphere-cone, = 10, a = 10°, 
<D = 90°. 



Figure 17.- Radial distribution of static pressure; 
x B /R n = 16.7, 15° sphere-cone, - 10, a = 10°. 


flow over an inclined body. 7 Therefore,, 
because of the entropy layer, the blunt- 
cone crossflow does not approach the 
pointed-cone distribution uniformly. 

The entropy layer is characterized 
by a nearly constant pressure, as is 
shown in figure 17. On the other hand, 
the density varies strongly on the wind- 
ward side of the body as may be seen in 
figure 18 for x^/R n = 16.7. The thick- 
ness of the entropy layer may be taken 
as the distance to the point where the 
density has a local maximum. This 
entropy layer is similar to that found 
in axisymmetric flow except that it 
develops faster (i.e., closer to the 
nose) on the windward side and more 
slowly on the lee side of the body. 

Smoothing, which was used to stabilize 
the calculation, causes the peak in density to be rounded off. A theoretical 
maximum for $ = 180° can be calculated by using the entropy corresponding to 
the minimum shock angle and by making use of the fact that the pressure is 
nearly constant. The theoretical maximum is about 15 percent higher than the 
calculated value. It is estimated that the smoothing had negligible effect on 
the density and crossflow angle for n/n greater than about 0.3 (see fig. 8). 



Figure 18.- 

V R n “ 

u - 10°, 


Radial distribution of density; 
16.7; 15° sphere-cone, M = 10, 


Comparison is made in reference 14 between inviscid and experimental 
(viscous) surface streamlines. 
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Comparison With Experiment 


Detailed experimental surveys of the flow field around blunted cones 
have been previously reported by Cleary (refs. 30 and 31). Experimental 
results from reference 30 are used to compare with the present theoretical 
results in order to confirm the validity of the proposed numerical methods. 
More complete comparisons between theory and experiment can be found in 
references 31 and 14, for both air and helium flows. 



Figure 19.- Surface pressure distribution; 15° 
sphere-cone, M = 10, a = 10°. 


Figure 19 presents the surface 
pressure distribution along each of the 
7 meridional planes along which the 
calculations were made. Theory and 
experiment are in excellent agreement 
although the experimental data tend to 
be slightly above the theory. This is 
consistent with the usual boundary- 
layer displacement effect. At 
x/R n = 20 the blunt-cone pressure is 
very nearly equal to the pointed-cone 
value shown on the far right of 
figure 19. 

Flow properties off the body 
surface are most easily studied 
experimentally by means of the impact 
or pitot pressure, as shown in fig- 
ure 20. The pitot-pressure distribution 
normal to the body is presented for all 
7 meridional planes and for an axial 
station of Xb/R n “ 16.7. The theory 
extends all the way to the shock in 
each case, while the data do not, except 
for $ = 60° and § = 150° where experi- 
mental shock positions are indicated by 
a sharp drop in pressure. 


In the viscous boundary layer the pitot pressure is low, going to zero at 
the body surface. Evidence of a very thick viscous layer is indicated by the 
experimental data for the lee side of the body, $ = 0 and 30°. On the remain- 
ing meridional planes the boundary layer is very thin. The entropy layer is 
clearly indicated in the experimental data by the peak in the pitot pressure 
just off the windward side of the body. The theory underpredicts the peak 
value because of the previously described smoothing employed in the calcula- 
tion. A theoretical maximum of Cp^ = 9.03 for $ = 180° is determined by the 

use of the total pressure at the calculated minimum shock angle (fig. 13). 
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Figure 20.- Shock-layer pitot-pressure distribution; 15° sphere-cone, M m = 10, 

a = 10°, x B /R n = 16.7. 


Aside from the noted differences due to the boundary layer and the entropy 
layer, the proposed numerical method appears to give an adequate prediction 
of shock- layer properties. 


CONCLUDING REMARKS 


Characteristics theory for three-dimensional steady flow was reviewed and 
the compatibility equations were derived in terms of pressure and stream 
angles as dependent variables. It was argued that the major practical dif- 
ficulty encountered in bicharacteristic methods results from the need for 
numerical differentiation and interpolation of randomly spaced data. Further- 
more, it was observed that shock- layer coordinates are essential to a simple 
treatment of circumferential derivatives on the shock surface, A reference 
plane method was therefore adopted in which an equal number of points are 
equally spaced between the body and shock along each reference plane. This 
mesh introduces the added complication of nonorthogonal coordinates into the 
equations but allows the use of simpler numerical techniques. 

With the constraint of a uniformly spaced mesh, particular characteristic 
lines are not followed as in more standard characteristics methods. The 
role of characteristics theory in the reference plane method is to determine 
how the finite difference equations are to be locally solved. Since the 
characteristics are not traced throughout the shock layer, the possible 
coalescence of waves to form embedded shocks is not determined during the 
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calculation. However, characteristic lines can be traced afterward, from the 
numerical solution. Neglecting embedded shocks is not a serious deficiency 
since the calculations are correct for weak shocks (entropy rise being third- 
order in the deflection angle) . In situations where strong shocks are sus- 
pected or known to occur (such as on flared bodies) special treatment would 
be required just as in standard two-dimensional methods. 

Results for the three-dimensional flow around inclined blunted and 
pointed cones were presented to demonstrate the methods described and devel- 
oped in this paper. Comparisons with a linearized characteristic method 
illustrated nonlinear angle-of-attack effects on crossflow parameters. 

Major effects of bluntness, which persist at large distances from the nose, 
were well predicted by the present methods. Predictions of surface and pitot 
pressures were found to be in good agreement with available experimental 
results for a 15° sphere-cone at 10° angle of attack. 

While present examples were limited to bodies of revolution, the methods 
are not so restricted, being generally applicable to smooth bodies without 
axial symmetry. More complicated shapes typical of high-speed aircraft will 
require additional development, primarily in the area of special boundary 
conditions. For example, it should be possible, with special treatment of the 
leading-edge boundary condition, to calculate the flow over wings with 
supersonic leading edges. Finally, while there are no inherent limitations in 
angle of attack, it is clear that the flow near the leeward meridian must 
become wake- like at a sufficiently large inclination. In this case, it is 
essential to include the effects of viscosity and possible secondary shocks. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, July 7, 1969 
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APPENDIX A 


DIRECTION COSINES FOR NONORTHOGONAL COORDINATES 


In the derivation of the compatibility equations for the reference plane 
method it was necessary to transform from streamline coordinates to a non- 
orthogonal system consisting of s*, n, and The direction s* is the 
projection of the streamline on the meridional plane, n runs from body to 
shock in the meridional plane, and ? is the out-of-plane direction (see 
fig. 3). This transformation is expressed as 

x. = e?.z.* (Al) 

i ID J 

where 

\ = (s,n,t) 

and 

Zj* = (s*,n,C) 

It is the purpose of this appendix to write out the expressions for the 
direction consines e| • . The analysis is made for the more general 
coordinates and is later specialized to the case where £ = s*. 

The directions of the nonorthogonal coordinates are determined 

step by step during the calculation as the shock shape is constructed (see 
section on Computational Procedure) . When the location of the shock is known 
in terms of x,r,$, the locations of field points between the shock and body 
are also determined by the prescribed spacing of the points along the body 
normals. The direction cosines can then be obtained by numerical differenti- 
ation as described below. This process is described for the shock wave but 
applies generally to each Z curve. 

Let the shock surface be given by 

r s = r s (x,$) (A2) 

The intersection of the coordinate surface g = constant with the shock 
surface defines a curved line in space (a £ curve; see fig. 3) . The x,r 
coordinates of this space curve depend on the meridional angle, and on 

the body station, x^ , which locates the £ = constant surface. For each body 
station this dependence is denoted by 

x s = x s (0) (A3a) 

r s = r s (« (A 3b) 


40 



Equations (A3a) and (A3b) define two angles 


tan Ax 


_ 1 _ 

r 

s 


dx 



d<I> 


(A4) 


1 ^ 

tan Ar = — (A5) 

s 



which can be evaluated numerically. 
(The Fourier method previously 
described is used.) In figure 21 the 
unit vector £ is related to e^ 
through two rotations, the first rota- 
tion by Ax and the second by 6. 

It may also be verified from this fig- 
ure that 


tan 6 


dr 
r d$ 


cos 


Ax 


and, by equation (A4) , 


Figure 21.- Coordinate vectors and angles. 


tan 6 = tan Ar cos Ax (A6) 


The remaining^two unit ^vectors are defined by angles X and a; n is rotated 
by X from e^ and £ is rotated by a from e x . Thus, one may write 

£ = cos a e + sin a e (A7) 

x r 

n = -sin X e + cos X e (A8) 

x r v J 

C = cos 6 sin Ax e + sin 6 e + cos 6 cos Ax e^ (A9) 

x r 0 v J 

The angle X, which determines the direction of n, may be arbitrarily 
specified. It is usually selected so that n is normal to the body surface. 
The angles a and 6 are determined at the shock by the calculated shock 
shape, and at the body by the given body shape. At intermediate points 
between the body and shock, a and 6 can be similarly calculated from the 
x,r,$ coordinates of the mesh points. These coordinates are known for any 
specified mesh spacing once the shock position is determined. 


It must be noted, however, that <5 is distinct from the angle 6 which 
is used in the shock conditions developed in appendix B. From figure 21 it 
is found that 


, _ dr - dx tan a 

tan 6 = - -- -jT’ — 

r d$ 

or 

tan 6 = tan Ar - tan a tan Ax 

When n is chosen to be normal to the body axis. Ax = 0 and 6 = 6. 


(A10) 
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Having specified the directions of £,n,S by equations (A7) , (A8) , and 
(A9) , the task of determining et . , appearing in equation (Al) , can now be 
completed. Equations (A7) , (A8)^and (A9) are more conveniently written with 
index notation as 




i 

ill 


where 


✓ cos a 

sin a 

0 


11 

*r-> 
1 ^ 

[ -sin A 

cos A 

0 



\cos 6 sin Ax 

sin 6 

cos 6 cos 


The inverse transformation may be written 



where 


X. = v. .z. 
1 13 3 


v. . = (v. 0 1 
ij ij 


(All) 


(A12) 


(A13) 


Although the matrix inversion can be performed numerically, for this problem 
it is more efficient to do the algebra beforehand, once^and for all. This is 
done by taking scalar products of equation (All) with Xj to obtain terms 
like z\ • Xx = \>n, z\ * X 2 = Vi 2 , and so on. On the other hand, the scalar 
products of equation (A13) with z_. give expressions such as 

z\ • Xi = v n + v 12 (? • n) + v 13 (? * ?)• 


Z'l " Xi = v^fri * t) + V12 + v 1 3 Cn • O 

and so on. In this way, one obtains nine equations for the nine unknowns 
. The solution of this set of equations by Cramer’s rule can be formally 
written as (see, e.g., ref. 18) 


v. . 
iJ 


v 7i f 7j 

D 


(A14) 


where f ^ . represents the cofactor matrix of the coefficients of . 
(, 


n 


V 


(n • c) 2 ] 

[d • 6) (n • i)-Cz • n)] 

[C5 

• n) (n • £)-(£ ■ 

• 0] 

f 21 

d-C€ ■ c) 2 ] 

[(? 

• n)(t * ?)-(n • 

■ 8)] 

fsi 

f 32 


[1-U • n) 2 ] 

/ 


(A15) 
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where 


f 





and the determinant D is given by 

D = 1 - (n • O 2 - (C • n) 2 - (5 • c) 2 - 2(i • n) (5 • cHn • S) (A16) 


Equation (A14) provides the direction cosines between g,n,? and 
x,r,$ coordinates. The corresponding relationships with s,n,t coordinates 
are obtained by use of equations (7), (8), and (9) which express the stream- 
line directions in terms of e^, e^, e^ and may be written 


where 


x. = t. .X. 
1 13 3 


cos <j> cos 0 cos $ sin 0 sin 


T ij =l 


-sin 0 


cos 0 


<-sin (j) cos 0 -sin <j> sin 0 cos 0 


The substitution of equation (A13) into (A17) yields 


(A17) 


(A18) 


x. 

l 


e . i z. 
lk k 


where c., is obtained from 
ik 


e.. - x . . v . , 
ik ij jk 


(A19) 


(A20) 


Equation (A20) gives the direction cosines between and s,n,t coordi- 

nates. The explicit expressions for these direction cosines are involved, so 
the final combination of terms indicated by equation _(A20) is left for the 
computer. The procedure is as follows. The matrix , defined by (A12), is 

calculated with equations (A4) , (A5) , and (A6) . The various scalar products 

in equations (A15) and (A16) are calculated from uj_ j , and the matrix v^j 

is then calculated from equation (A14) . Finally, is determined from 

equations (A18) and (A20) . 

The procedure described applies to general nonorthogonal coordinates 
The transformation used in equation (36) is a special case where the 
E, direction is the streamline projection s*. In this case, equation (A9) is 
replaced by 
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(A21) 


s * = cos 0 e + sin 0 e 


Then the subsequent relations will be modified accordingly and there is 
obtained 


e 


* 

ik 


T. .vt, 

ij jk 


(A22) 


which determines the desired transformation relationship indicated by 
equation (Al) . 
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APPENDIX B 


SHOCK -BOUNDARY CONDITIONS FOR THREE-DIMENSIONAL FLOW 


Consider an elemental portion of the shock surface with an outer unit 
normal N, as shown in figure 22 , and with unit vectors T and L completing 
an orthogonal set. The tangent vector T can be chosen such that the N-T 

plane is parallel to the direction of the 
free-stream velocity. This choice will 
permit evaluation of the jump conditions 
in the N-T plane with two-dimensional 
shock relations. Let s be a unit vec- 

oo 

tor parallel to the free-stream velocity 
vector and with components 

A A /V 

s = cos a e Y + sin a cos $ e_ 

oo A -L 

- sin a sin $ (Bl) 

The desired tangent vector T can be 
constructed from s and N by means of two vector- or cross-products. The 
first product 

aL = s x N (B2) 

oo 

produces a vector parallel to L, and the second product 

T = NxL (B3) 

results in a vector which is normal to both N and L, and which lies in the 
s^-N plane. The factor a in equation (B2) is equal to the sine of the 
angle between s^ and N and may be evaluated from the scalar product 

b = s * N = cos (s ,N) (B4) 

to give 


a = -s/l - b 2 = J 1 - (s M ■ N) 2 (B5) 



Figure 22.- Shock-wave normal and tangent 
vectors . 

in a cylindrical coordinate frame. 


Combining equations (B2) and (B3) and expanding the resulting vector triple 
product one obtains 


or 


T = I [Nx (s^ x N)] ' 
T = t - bN) ‘ 


(B6) 
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In terms of components T x , T^, and T$, equation (B6) becomes 


T 



(B7) 


Equation (B7) permits evaluation of the true inclination of the shock surface, 
and therefore allows the jump conditions to be determined with standard planar 
shock relations (see, e.g., ref. 32). The angle a* between and T is 
given by 


cos o f = s • T 

CO 


or 


cos a 1 


s T 

X X 


+ s T + s T, 
r r $ $ 


(B8) 


In the following development a* is expressed in terms of two angles 
measured in the cylindrical coordinates. 


Let o be the angle between e x and the trace of the shock surface on 
the plane $ = constant, and let 6 be the angle between e$ and the shock 

trace on the- plane x = constant as illus- 
trated in figure 23. The shock angle 6 is 
obtained by numerical differentiation as 
described in appendix A (see eq. (A10) and 
fig. 21). It is easily verified that the fol- 
lowing relations hold between the angles 
shown in figure 23. 




Shock 



<t> = Constant 


N = 
x 


-N sin a 
m 


N = N cos a 
r m 


N, = -N cos a tan 6 
$ m 


(B9) 

(BIO) 

(Bll) 


where N m is the projection of N on the 
x-r plane and satisfies the relation 


Figure 23.- Shock angles. 


N 


m 


= 

$ 


Substitution of from equation (Bll) gives 

N = — 1 - (B12) 

m Jl + cos 2 a tan 2 6 
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and the shock normal vector may finally be written in terms of a and 6 as 

(-sin a e + cos a e - cos a tan 6 e_) 

a Y Y* ffl 

N = X . - — (B13) 


1 + cos 2 a tan 2 6 


The true shock angle can now be evaluated in terms of a and 6 by equations 
(Bl) , (B7), (B8) , and (B13) , and the jump conditions can be calculated. 


It is now necessary to determine the flow angles 0 and $ measured 
relative to the meridional planes. 



The streamline direction, measured in 
the N-T plane (fig. 24), is turned from the 
free stream by the angle 02*. Since the 
streamline tangent §2 lies, by definition, in 
the N-T plane, one may write 

s 2 = -sin(a l - 0 2 f ) N + cosOr* - 0 2 ! )T 


Figure 24.- Shock and flow angles in the 
n-t plane. 


(B14) 


Using equations (B7) and (B13) , the vectors N and T can be written in terms 
of their components to give 


s 2 = (AN x + BT^ + (AN r + BT> r + (AN $ + BT $ )S $ (BIS) 

where 

A = -sin (a T - ©2 * ) 

B = cos (a f - 02 1 ) 

Equation (7), on the other hand, gives S 2 in terms of 0 and <J> as 
follows : 

s 2 = cos <j >2 cos ©2 e x + cos (j > 2 sin 0 2 + sin 2 e § (B16) 

Thus, by equating components of equations (B15) and (B16) one obtains finally 


and 


tan 0 


-sin (a 1 - © 2 1 + cos (a 1 - 02 r )T^ 

-sin(a f - 02 ! )N x + cos (a * - 02 f )T 


(B17) 


sin <j> = -sin(a f - 02 f )N^ + 

If for given free-stream conditions, p^, p^, 
conditions are written in the form (see ref. 


cos (cr 1 - 0 2 f )T o (B18) 

V^, the standard planar shock 


32) : 
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P = pO') 
P = p(o') 

e = 0 (o 1 ) 


(B19) 


then the equations developed in this appendix allow the three-dimensional 
shock conditions to be functionally written as 


p = p(a;6,$,a) 
p = p(a;8,<I>,a) 
0 = 0 (a ;<$,$, a) 

<f> = <t>(o '>&>$, °0 


(B20) 


This is the form of the shock conditions employed in equation (62]. 

The overall shock calculation is performed in a straightforward iterative 
manner. This procedure is started with known field data, including a and 5, 
on an initial data surface. A shock point on the new data surface is 
determined by the average value of a between the initial and new shock 
points , 

a = 1/2 (a x + a 2 ) 

To begin, a 2 is set equal to o\ and then subsequent estimates are made 
(either by the Newton method or by the bisecting method) until the pressure 
from equation (B20) agrees with the pressure calculated from the compatibility 
equation (57) . During each iteration the current value of a = cr 2 and the 
old value of 6=6} are used to determine the true shock angle a T by 
equation (B8) . The pressure and other variables are then calculated from 
equations (B19) . When this has been done for all the shock points on the 
new data surface, new values for 6 can be determined by numerical differen- 
tiation with respect to $ as described in appendix A. The entire process 
can then be repeated to refine the results; this is discussed in the section 
on Global Iteration . 
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APPENDIX C 


SURFACE BOUNDARY CONDITIONS FOR BODIES WITHOUT AXIAL SYMMETRY 


For noncircular bodies the surface boundary is complicated by the fact 
that the surface normal is not in the meridional plane. This means the 
boundary condition will involve both flow angles 0 and <f>. (In terms of 
velocities, all three components, u, v, w, enter into the conditions since 
none of the components are parallel to the body surface.) In this appendix 
the boundary condition for 0, which was given previously in equation (61a), 
is derived from the tangency condition on the velocity vector. 

Let the equation of the body be given by 

g(x,r,<f>) = r - f(x,0) (Cl) 

The unit outer normal to the surface may be expressed 

N = (C2) 

I Vg | 


where V is the vector gradient operator. In terms of cylindrical 
coordinates, one obtains 


9f - - 1 9f a 

3x 6 x + 6 r r 30 6 0 


1 + 


(MY + (I!£Y 

\ 3x / \r 30/ 


(C3) 



The derivatives of f in equation (C3) 
are related to the surface inclination 
angles, figure 25, according to 


' r A 

^ \\ 


tan 6 


{<£ constant) 

Figure 25.- Surface inclination angles. 


and 


tan <5 


0 


■ (S). 

• i ffi) 


(C4) 


(C5) 


From equation (7) the streamline direction is expressed in terms of 
flow angles 0 and <f> 


s = cos <b cos 0 e + cos d> sin 0 e + sin <f> e. 

x r 4> 


CC6) 
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The tangency condition, s * N = 0, may now be obtained from the scalar product 
of equations (C3) and (C6) . This gives 

cos <j) cos 0 + cos <j> sin 0 + sin <|> = 0 (C7) 

where N x , N r , are the components of N along e x , e r , e$, and are 

identified by equation (C3) . Equation (C7) may be written as a quadratic 
relation in tan 0 

(N^ 2 - N^ 2 tan 2 <j>)tan 2 0 + 2 ^ x ^ r tan e + C^ x 2 “ N^ 2 tan 2 <j>) = 0 (C8) 

The solution is 


tan 0 



[1 - (N $ /N r ) 2 tan 2 <p] 


(C9) 


Rewriting equation (C9) in terms of the surface inclination angles given 
by equations (C4) and (C5) , one obtains 


tan 0 = 


^tan + tan 6^ tan cj> J 1 + tan 2 6^ - tan 2 6^ tan 2 c()^ 


(CIO) 


(1 - tan 2 6^ tan 2 <p) 


where the positive root is chosen so that 0 is decreased when <j> < 0 and 

6 * > 0 . 

$ 

Equations (CIO), (C4) , and (C5) determine the flow angle 0 in terms of 
the crossflow angle $ and the body geometry. It is easily verified that 
for zero crossflow, cf> = 0, equation (CIO) reduces to 

tan 0 = tan 6 

x 


which is the usual condition for circular bodies. 
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TABLE I.- EFFECT OF GLOBAL ITERATION ON SHOCK ANGLE 
[15° sphere-cone; M = 10, a = 10°, 15 points and 7 planes] 

(a) Calculation from x/R n = 2.0 to x/R n = 3.0 


P 1 ane , 

Iterations 

0 

i 

deg 

a, deg at 

x/R n = 3.0 

0 

29.0763 

29.0763 

30 

27.6271 

27.6275 

60 

23. 7124 

23.7137 

90 

18.8404 

18.8417 

120 

15.5193 

15.5194 

150 

14.2633 

14.2633 

180 

13.9813 

13.9813 


(b) Calculation from x/R n = 10.0 to x/R n = 11.0 


P 1 ane , 

Iterations 1 

0 

1 

deg 

a, deg at 

x/R n =11.0 

0 

21.2083 

21.2083 

30 

19.5249 

19.5273 

60 

16.1617 

16.1651 

90 

16.0245 

16.0260 

120 

17.4741 

17.4756 

150 

17.4404 

17.4407 

180 

17.4416 

17.4415 



St 




TABLE II.- ACCURACY AND COMPUTING TIME 
[15° sphere-cone; M = 10, a = 10°. Calculation from x/R n = 2 to x/R n = 3.] 


(a) Shock angles and surface pressure on $ = 90° plane, x/R n = 3 







Planes 







K = 3 



K = 5 


K = 7 

Points 

0 

6 

Pb 

a 

6 

Pb 

a 

6 

Pb 

J = 5 

19.3114 

-10.3313 

840.42 

19.5919 

-10.1783 

834.02 

19.2893 

-10.3555 

828.85 

J = 10 

18.9847 

-10.3961 

861.15 

18.9337 

-10.5549 

849.52 

18.9318 

-10.5644 

849.36 

J = 15 

18.9138 

-10.4090 

864.56 

18.8482 

-10.6017 

853.39 

18.8404 

-10.6294 

853.44 


(b) Computing time 







Planes 






K = 

3 


K = 

5 


K = 

7 

Points 

Total 
points , 
N 

Time, 

t, 

min 

Unit 

time, t/N 

Total 
points , 
N 

Time, 

t, 

min 

Unit 

time, t/N 

Total 
points , 
N 

Time , 
t, 

min 

Unit 

time, t/N 

LO 

n 

90 

0.36 

0.400xl0 -2 

150 

0.57 

0.380xl0 -2 

210 

0.7 

0.333x10-2 

J = 10 

420 

.75 

.179xl0 -2 

700 

1.17 

. 167xl0- 2 

980 

1.70 

. 174x10-2 

J = 15 

1035 

1.35 

.130xl0- 2 

1725 

2.30 

.133x10-2 

2415 

3.34 

.138xl0- 2 
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TABLE III.- 15° POINTED CONE SOLUTION 


PERFECT GAS 

GAS CONSTANT = 0.17160E 0* GAMMA * O.IAOOOE 01 
FREE-STREAM CONDITIONS 

H= Q. 106006 02 V= 0-396626 04 P* O.IOOOOE 01 RHO* 0.100006-04 T * 0.5S275E 02 H* 0.360006 06 HT=* 0.821526 07 

EOUISPACEO STARTING DATA 11 POINTS 9 PLANES 

PLANES EQUALLY SPACEO 

STARTING DATA NORMAL TO BODY SURFACE 

ANGLE OF ATTACK * 10.00 DEG 

CHARACTERISTIC METHOD 

I N( 21 ) = 0 COR» N22 « 0 ITR, INI25)=i 2 BC , INXI61- 3 BODY DENSITY 

FL f A ) = 0.1800E 01 STEP SIZE, FL(0)= O.IOOOE 01 SMOOTH 


PLANE 1 ANGL = 0.00 DEG 

LEEkARO PLANE 

FIELD DATA 

SHOCK ANGLES, DEG SIGMA= 18.5LZ4 DEL T A= -0.0000 

X R TH6TA PHI P RHO H V M M* 

0 * 983956 DO D.32785E 00 0.24764E 00 -O.OOOOOE-38 0.27055E 01 0.19796E-04 O.47036E 06 0.39337E 06 0.89927E 01 0.89927E 01 


0.9 8555k 00 O’. 32 1 8 6E 00 0.2494SE 00 -O.OOOOOE-38 0.27456E 01 

0.98716E OO 0.31587E 00 0.25107E 00 -O.OOOOOE-38 0.27781E 01 

0.98076(6 00 0.30988E 00 0.25257E 00 -O.OOOOOE-38 0.280476 01 

0.990376 00 O.30389E 00 0.25400E 00 -O.OOOOOE-38 0.28266E 01 

Q.99197E 00 D.29790E 00 0.25543E 00 -O.OOOOOE-38 0.284456 01 

0.993586 00 0.29191E 00 0.256906 00 -O.OOOOOE-38 0.28589E OL 

0.99516E 00 O.20592E 00 0.25844E 00 -O.OOOOOE-38 0.28701E 01 

0.996796 00 0.27993E 00 0.259946 00 -O.OOOOOE-38 0.28786E 01 

0.99839* 00 0.27394E 00 0.26LL6E 00 -O.OOOOOE-38 0.288556 01 

O.IOOOOE 01 0.26795E 00 0.26180E 00 -O.OOOOOE-38 0.28922E 01 

PLANE 2 ANGL * 22.50 DEG 

FIELD DATA 

SHOCK ANGLES, DEG SIGMA= 18.5966 DELTA- 1.3917 

X R THETA PHI P 

0.98339E 00 0.32995E 00 0.2489BE 00 -0.69236E-01 0.32071E 01 

0.98505E 00 0.32375E 00 0.24930E 00 -0.69743E-01 0.3176TE 01 

0.98671E DO 0.31755E 00 0.24990E 00 -0.70020E-01 0.31529E 01 

0.988376 00 0.31U5E 00 0.25075E 00 -0. 700216-31 0.31339E 01 

0.99033E 00 0.30515E 00 0.25182E 00 -0.69673E-01 0.311B2E 01 

0.99169E 00 0.29895E 00 0.25312E 00 -0.6B884E-Q1 0.31043E 01 

0 • 993356 00 0.29275E 00 0.25461E 00 -0.67720E-01 0.309116 01 

0.99501E 00 0.28655E 00 0.25629E 00 -0.66943E-01 0.30774E 01 

0.99668E 00 0.280356 00 0.258IIE 00 -0.68538E-01 0.30619E 01 

0.99834E 00 0.27415E 00 0.26000E 00 -0.74793E-D1 0.30438E Oi 

O.IOOOOE 01 0.26795E 00 0.26180E 00 -D.85422E-01 0.30238E 01 


PLANE 3 ANGL = 45.00 DEG 

field DATA 

SHOCK ANGLES, OEG SIGMA- 18.5775 DELTA- -1.5278 

X R THETA PHI P 

0 .983 3 IE 00 0.330256 00 0.255046 00 -0.12297E DO 0.50977E 01 

0.98497E OO 0.32402E 00 0.25413E 00 -0.12138E 00 0.49424E 01 

0.98664E 00 0.31779E 00 0.25357E 00 -0.12055E 00 0.48026E 01 

0.988316 00 0.31156E 00 0.25335E 00 -0.11942E 00 0.46769E 01 

0 . 98998 E 00 0.30533E 00 0.253496 00 -0.11793E 00 0.45640E 01 

0.99165E 00 0, 2991 OE 00 0.253996 00 -0.11604E 00 0.64624E 01 

0.99332E 00 D,29287E 00 0.254B5E 00 -0.11385E 00 0.437086 01 

0.99499b 00 0.28664E 00 0.256D8E 00 -0.11203E OO 0.42875E 01 

0 » 99666E 00 0.28041E 00 0.2S766E 00 -0.H23TE 00 0.42097E 01 

0.99633E 00 O.27410E 00 0.25959E 00 -0.11741E 00 0.41346E 01 

O.IOOOOE 01 0.267956 00 0.26280E 00 -0. 12786E 00 0.4D623E 01 

PLANE 4 ANGL » 67.50 DEG 

FIELO DATA 

SHOCK ANGLES, DEG SIGMA* 28.2994 DELTA- -2.5439 

X R THETA PHI P 

0.98464E OO 0 * 32 52 66 00 0.25L40E 00 -0.158156 00 0.82016E OL 

0.98618E 00 0.31953E 00 0.25099E 00 -0.15560E 00 0.79933E 01 

0.98771E 00 0.3138 OE OO 0.25086E 00 -0.15280E 00 0.78012E 01 

0.989256 00 0.30806E 00 0.25103E 00 -0.14971E 00 0.76251E 01 

0.99079E 00 0.30233E 00 0.251506 00 -0.14630E 00 0.74638E 01 

0.99232E 00 0.29660E 00 0.25228E 00 -0.14252E 00 0.73160k 01 

0 . 993866 OO 0.290876 00 0.25339E 00 -O.13045E 00 0.71800E 01 

0.99539E 00 0.28514E 00 0.254846 00 -0.13460E 00 0.70548E 01 

0.99693E 00 0.27941E 00 0.256706 00 -0.132306 00 0.69389E 01 

0.998456 00 0.27368E 00 O.25901E 00 -0.133106 00 0.68305E OL 

0.10000* 01 D.26795E 00 0.26180E 00 -0.137056 00 0.67275E 01 


0. 19989E-04 0.4B074E 06 0.39331E 04 0.896906 01 0.89690E 01 
0. 20145E-04 0.48267E 06 0.393266 04 0.89499E 01 D.89499E 01 
0.20269E-04 O.40432E 06 0.39321E 04 0.89337E 01 0.89337E 01 
0.20366E-04 0.48577E 06 0.39318E 04 0.89196E 01 0.89196E 01 
0 • 20441E-04 0.48704E 06 0.39315E 04 0.89072E 01 0.89072E 01 
0.205056-04 0.48799E 06 0.39312E 04 0.88980E 01 0.889806 01 
0. 205706-04 0.48836E 06 0.39311E 04 0.8B944E 01 0.88944E 01 
0 . 20640E-04 0.48813E 06 0.39312E 04 0.889666 01 0.889666 01 
0.207 06E-04 0.487746 06 0.39313E 04 0.890046 01 0. 890046 01 
0.207626-04 0.48757E 06 0.393136 04 0.890216 01 D. 890216 01 


RHO H V M M* 

0,219866-04 0.51055E 06 0.39255E 04 0.86864E 01 0.86659E 01 

0.21793E-04 0.51019E 06 0.39256E 04 0.868976 01 0.86689E 01 

0.21615E-Q4 0.51052E 06 0.392556 04 0.86867E 01 0.86657E 01 

0.21439E-04 0.511636 06 0.39252E 04 0.86767E 01 0.86557E 01 

0.212396-04 0.51385E 06 0.39246E 04 0.86567E 01 0.863606 01 

0.20962E-04 0.51832E 06 0.39235E 04 0.861686 01 0.859666 01 

0.204476-04 0.52910E 06 0.39207E 04 O.0S225E 01 0.850326 01 

0.193096-04 Q.557B1E 06 0.39134E 04 0,828486 01 0.82665E 01 

0.17025E-04 0.62948E 06 0.389516 04 0,776236 01 0,774446 01 

0. L35 796-04 0.78452E 06 0.38550E 04 0.688176 01 0.686296 01 

0. 102026-04 0.10374E 07 0.37889E 04 0.588176 01 0.586096 01 


RHO H V M M* 

0.2B462E-04 0.626B7E 06 0.389576 04 0.777986 01 0.772306 01 

0.27594E-04 0.62689E 06 0.389576 04 0.777976 01 0.772346 01 

0.267396-04 0*628646 06 0.38953E 04 0.776806 01 0.771256 01 

0.258796-04 0.63253E 06 0.389436 04 0.774206 01 0.768786 01 

0.249916-04 0.639186 06 D.38926E 04 0.76983E 01 0.764576 01 

0.240436-04 0.64960E 06 0.388996 04 0.763116 OL 0.758066 01 

0. 229856-04 0.66556E 06 0,308586 04 0.753116 01 0.74831E 01 

0.21705E-04 0.69138E 06 0.387916 04 0.73764E 01 0.733106 31 

0.199136-04 0.73992E 06 D.38666E 04 0.71073E 01 0.706336 01 

0.17061E-04 0.848186 06 0.38385E 04 D.65900E 01 0.654566 01 

D.12596E-04 0.11287k 07 0.376476 04 0.56028E 01 0.55585E 01 


RHO H V M M* 

0.353556-04 0.81193E 06 0.38479E 04 0.675216 01 0.66697E 01 

0.34195E-04 0.81815E 06 0.38*636 04 0.67236E 01 0.664406 01 

0.33049E-04 0.826176 06 0.3B442E 04 0.668726 01 0.66109E 01 

0.31909E-04 0.8 363 7E 06 0.384166 04 0.66417E 01 0.65690E 01 

0 . 30759E-04 0.84930E 06 0.383826 04 0.65852E 01 0.65164E 01 

0.295666-04 0.86607E 06 0.383386 04 0.651376 01 0.64491E 01 

0.28239E-04 O.08991E 06 0.38276E 04 0.641546 01 0.635546 01 

0 . 26552E -04 0.92994E 06 0.381716 04 0.625866 01 0.62D34E 01 

0.241 75E-04 0.10046E 07 0.37975E 04 0.599076 01 0.59398E 01 

0.21092E-04 0.11335E 07 0.376346 04 0.55892E 01 0.55413E 01 

0.18061E-04 0.130376 07 0.371796 04 0.514856 01 0.510206 01 
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TABLE III.- 15° POINTED CONE SOLUTION - Concluded 



PLANE 5 ANGL = 90.03 DEG 

FIELD DATA 

SHOCK ANuLEStr- UEG SIGMA= lb. 01^6 DELTA* -2.3858 

X R THETA PHI P 

0 . 98591 t 03 0.32052E 00 0.24054C 00 -0.173246 GO Q.12522E 02 

0.98732E 00- 0.31526E 00 0.24160E 00 -0. 169I5E 30 0.I2337E 02 

0.988736 3t> 0.3I000E 00 0.242876 00 -0.16489E JO 0.12162E 02 

0 • 7931 46 00 0.30475E OO 0.244336 00 -0.16Q42E 30 0.U996& 02 

0.99I55E 00 3.299496 00 0.24599E 00 -0.15575E GO O.U038E 02 

0 . 99296.fr 00 0.29423E 00 0.247866 00 -0.15G85E GO 0. 11688E 02 

0 .994366 00 0.28898E 00 D.24997E 00 -D. 14568E 00 0.11544E 02 

0.995778 00 0.28372E 00 0.25235E 00 -0.14326E 00 0.U408E 02 

0.997186 GO Q.27846E 00 Q.25506E 00 -0.135036 00 0.11277E 02 

0.9985*6 90 0.27321E 30 0.258I8E 00 -0.131216 00 0.111516 02 

0. 100006 01 0.267956 30 0.261806 00 -0.129656 GO 0-U330E 02 

PLANE 6 ANGL *112.50 DEG 

FIELD DATA 

SHOCK ANGLES. DEG SIGMA* 17.7465 OfcLTA* -1.9369 

X R THETA PHI P 

3.987196 33 0.31577E 30 0.226J4E 00 -9.16273E 03 D-17628E 02 

0.988476 GO 3.31099E 00 0.22929E 03 -0.15826E 00 0.17593E 02 

0.989756 00 0.306206 00 0.23230E 00 -0.153656 00 0.175476 02 

0.991036 00 0.30142E 0 0 0.23S39E 03 -0.148866 03 0.17492E 02 

3.992316 30 0.29664E 00 0.23860E 00 -3.143866 00 3-174316 02 

3.99359E 30 0.2918&E 00 0.24193E 03 -0.13859E JO 3.173636 02 

0.99487E 30 0.287086 00 0.24S43E 00 -Q.13299E 00 0.17288E 02 

J. 99616-6 00 0.28229E 00 0.24913E 00 -0. 12706E 00 0.172I0E 02 

0 • 997446 00 0.277516 00 0.2533SE 00 -0.123956 00 3.17L27E 02 

0.99872b 30 0.272736 00 0.257246 00 -0.11478E 00 0.17338E 02 

0.100306 Oi 0.26795E GO 0.26180L 33 -0.108456 GO 3.16937E 02 

PLANE 7 ANGL =135.00 DEG 

6 i £L0 DATA 

SHUCK ANGLES. Dfc G SIGMA* 17.5849 DELTA* -1.1334 

X R THETA PHI P 

J. 987926 GO 0.3L302E 00 0.21284E 03 -0.12801E 00 0.226086 02 

0.989136 00 3.308526 00 0.217586 30 -0.124216 JO 3.22765E 02 

0.990346 00 0.304016 00 0.22229E 03 -0.120276 uO 0.22893E 02 

D. 991546 00 3.299506 30 0.22700E 03 -0.116176 00 3.229956 32 

0.99275 L 00 0.29499E 00 0.23172E 00 -0.111876 DO D. 233736 02 

0.993966 03 0.29049C 00 0.236496 00 -0.107356 30 0,231296 02 

0.995176 30 0.285986 00 0.241326 00 -0.102536 00 3.23165E 0? 

0.996386 30 D.28147E 30 0.24623E 03 -0. 972726 -01 0.23179E 02 

0.99758E 30 0.276966 OO 0.251266 00 -0.91283E-01 0.231726 02 

0.998796 30 0.272466 OO 0.256446 00 -0.845106-01 0.231446 32 

0.130306 01 0 . 26795L 00 0.261806 00 -0.773806-01 3.23097E 02 


PLANE 8 ANGL =157.50 DEG 

6 1ELD DATA 

SHUCK ANGL £ S ► DEG SIGMA* 17.4726 DELTA* -3.6073 

X R THETA PHI P 

0.988456 00 0.311056 00 0.20345E 00 -3.700286-01 3.263036 02 

0.98963'E GO 0.306746 00 0-.20951E 03 -0.6785TE-31 3, 266416 32 

0 . 993 76t GO 0.302436 30 0.216456 03 -0.6560&6-J1 0.269356 02 

0.99L91t 00 0 . 290 1 2fc 00 0.221316 03 -0.632676-01 0.271896 02 

0.993376 00 0 . 29 38 1 E OO 0.22712E 03 -3. 608276-01 0.274046 02 

3.994226 30 0.289506 30 0.232906 00 -0.582636-01 0.275846 02 

0.995386 00 0.265196 00 0.236b8fc 30 -0.555276-01 0.27728E 02 

0 . 996 5 3 fc 00 0.280886 00 0.24445E 33 -0. 525416-01 0.27837E 02 

0.997696 30 3.276576 00 0.250246 03 -0.491966-01 3.279126 02 

0.998846 00 0.272266 00 0.25601E 00 -0.453136-01 3.27955E 02 

0.100006 Gl 0.26795E 00 0. 261806 OG -3.408846-01 3.279596 02 

PLANE 9 ANGL =180.00 DEG 

6 tf i. 0 DATA 

S*. -K ANGLES, DEG SIGMA* 17.4467 DELTA* 0.3000 

x ft th&ta phi p 

0.98854E 00 0.31070E 00 0.203196 00 -0.587506-09 0.27683E 02 

0.98 96 96 00 0.306436 OO O'. 206666 00 -0.571746-09 0.280876 02 

0.993836 00 0.30215E 00 0.213006 00 -0.554856-09 3.28444E 02 

0.991986 00 0.297886 00 0.21923E 03 -0.53696E-J9 0.287566 02 

0.993136 00 0.293636 00 0.225386 03 -0.51904E-J9 3.290256 02 

J.99427E 37 0.289336 00 0.23149E 00 -0.502346-09 0.292526 02 

D. 995426 30* 0.28505E 00 0.237566 00 -0.487576-09 0.29439E 02 

0.99656€ 00 0.2B078E 00 0.24361E 00 -0.47286E-09 0.295876 02 

0.99771E OO 0.27650E 00 0.24967E 00 -0. 451226-09 0.29694E 02 

0.998866 00 0.27222E 00 0.25574E 00 -0.4L598E-09 0.29762E 02 

0. 1 OOOOt 01 0.267956 00 0.26180E 00 -0.37021E-09 3.29793E 02 

THETA BODY *0.2617996 00 


RHO H V M M* 

0.411346-34 0.1D663E 07 0.37812E 34 G.?7899fc 31 3.57*586 31 

0.40064E-04 0.107786 07 3.37782E 04 0.57542E 01 3.567446 31 

0. 39019E-D4 3.1091GE 07 0.37747E 04 3.571416 01 0.5638BE 31 

0.37958E-04 0.H061E 07 0.37707E 04 0.56687E 01 0.5598tE 31 

0 . 3 1> 8 686-04 0.U238E 07 0.37660E 04 3.56L69E 01 0.5551‘E 51 

0.35739E-04 3.11446E 07 9.37605E 04 0.55576E 01 0.54964E 31 

0 . 34566 E -0 4 0. 1L689E 07 3.37540E 04 0. 54899E Ol 3.543366 31 

0.3331 7E-04 0.U984E 07 0.37461E 04 0.54137E 01 3.53593E 31 

0.31777E-34 9.1 242 1 E 07 0.373456 04 0.529816 Ol 0.52515E 31 

0.293916-34 0.13279k 07 0.37114E 04 3.5:9246 91 0.5D532E 3L 

0.2571 IE-04 0.150156 07 0.36643E 04 0.472826 01 0.46933E 31 


RHO H V M m* 

9.451876-34 G.13654E 07 0.373136 34 0. 503836 01 3.4944R6 31 

0.44637E-04 0.13795E 07 0.369756 04 3.497766 01 0.49177E 01 

0.44049E-34 0.13942E 07 0.36935E 04 0.49459E Ol 3.488986 31 

0 . 434286-04 0.14097E 07 3.368936 04 0.491306 Ol 0.486076 31 

0.42774E-04 0.14263E 07 O.36048E 04 0.48784E Ol 0.483306 31 

0.420766-04 0.14443E 07 0.36799E 04 0.484L5E OL 0.47973E 01 

0.4 1 296E -0 4 0. 14653E 07 3.36742E 04 3.47993E Ol 3.47586E 31 

0.40302E-04 0.14946E 07 0.366626 04 0.474166 31 3.473536 31 

0.388576-04 0.154276 37 3.36531E 04 0.465036 Cl 0.46179E 31 

0.36921E-04 0.161516 07 0.363326 04 0.45202E Ol 0.44919E 31 

0.34927E-34 3.16973E 37 0.36lQ5fc 04 3.438196 5 1 0.435156 31 


RHO H V M M* 

0.477666-04 0.165666 07 0.362186 04 0.44492E 31 3.4414 76 j1 

0.47737E-34 0. 166916 37 3. 361836 34 3.4428 36 .31 3.439586 31 

0.476526-04 0.16814E 37 3.36149E 04 3.443786 Ol 3.43775t 31 

0. 475126-04 3.16939E 07 9.361146 34 C.43874fc Gl 0.4359^6 31 

0.473166-04 0.170676 07 0.36379E 04 3.436666 01 3.534966 31 

0.47062E-34 0. 1 72316 37 ''. 363426 04 3.43451E 01 3.432 L4E 31 

0.467516-04 J.17342E 07 0.363036 04 0.432276 01 D.4331IE 31 

0.463876-04 0.174896 Jl 0.359626 34 3.429966 01 3.42d93t 31 

0. 459096-04 0.176666 07 3.359136 34 0.427226 01 3.425646 51 

0.45049E-04 0.179816 07 0.358256 34 0.422426 01 0.421396 31 

0.435896-34 0.185466 37 9.356676 04 0.414116 01 3.41274E 31 


kHO H V M M* 

0. 491656-34 0.187256 07 3.356176 04 0.411546 31 3.41059t 31 

0.495326-34 0.168256 37 0.355896 04 3.413126 01 0.4j923t 91 

0.498376-34 0.189166 07 0.355636 04 0.408846 01 0.438316 01 

3.503836-34 0* 190331 37 3.355396 04 0.407666 31 3.436896 31 

0.502 726 -04 0.190796 0 7 3. 3551 76 04 3.436566 Ol 9.43->85& 31 

0.504046-04 3.191546 07 0.354966 04 0.405536 Ol 3.434886 31 

0.50474E-D4 0.192276 07 0.354756 04 0.404526 01 3.433946 31 

0.534626-04 0.193376 C7 3.354536 34 0.403426 31 0.4029'c 31 

0.503346-04 0.194396 07 3.354246 04 0.4G2346 31 3.431586 31 

0.501296-34 0.195186 07 0.353936 04 0.403566 01 3.4J3186 31 

0.499636-04 0.195866 07 0.353746 04 3.399656 Cl 3.399346 31 


RHO H V M M* 

0 . 496096 -04 0.195316 07 3.353906 C4 3.403396 01 3.403196 31 
0. 501146-04 0.196176 07 0.353656 04 0.399246 01 3.399246 31 
0.50556E-04 0 . L 9692 E C7 0. 353446 04 0.39824E 01 0.398246 31 
0.539406-04 0.19758E 97 0.353256 04 0.39737E OL 3.397376 Jl 
0.512966-04 0.198166 C7 0.353096 >4 0.396606 01 3.796676 31 
0.515356-04 0.198676 07 9. 35295k 04 0.395936 Oi 3.39593c ~\ 
0.51743E-04 0.199136 07 3.35281E 04 0.395316 Cl 3.395316 31 
0.518866-04 3.199S8E 07 0.35269E 04 0.394736 01 3.394736 01 
0.519836-04 0.199936 07 0.35269E 04 0.394276 OL 0.39427E 01 
0. 521006-04 0.199946 07 0.35259k U4 3.394266 01 3.394266 31 
0.52282E-34 0.19945E 07 3.35272k 04 5.3949CE 01 0.394936 31 
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